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I.  Synopsis 


As  recommended  in  the  brochure,  "Administration  of  U.S.  Air  Force 
Grants  and  Cooperative  Agreements  for  Basic  Research,"  this  report 
provides,  in  a  single  document,  a  permanent  record  of  the  progress  and 
significant  accomplishments  achieved  in  the  performance  of  the  research 
effort.  This  first  section  contains  summaries  of  research  background, 
research  objectives,  status  of  research,  publications  and  professional 
personnel  associated  with  the  research  effort.  Sections  II,  III,  IV 
and  V  are  technical  papers  resulting  from  the  research  effort  and 
Section  VI  is  a  Master  of  Science  thesis,  supervised  by  the  principal 
investigator,  which  contains  mathematical  details  and  several  computer 
codes  related  to  the  four  technical  papers  contained  in  this  report. 


»( 
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A.  Background 

The  research  results  reported  here  are  direct  extensions  of  earlier 
research  conducted  by  the  principal  investigator  under  the  auspices  of 
the  1979  USAF-SCEEE  Summer  Faculty  Research  Program  sponsored  by  AFOSR 
contract  No.  F49620-79-0038.  The  final  report  summarizing  the 

1/  *" 

summer  research  effort,  dated  17  August  1979  and  authored  by  the 


principal  investigatorV'itfhatttitl^i?  "A8£pt1v!*$i§S§i’ Processing  for 

Vv  •  i<V.  v<  ' 

Array  Antennas."  That^ ceoearah  Ua&Apehf§4ttiei:  St'fchd  Electronics 
,  ./  t  m*  n 

Systems  Division/Hanscom  Air  Forc^-BTa#e  :fiM  tf^oITabbration  with 


USAF  Research  Colleague,  Dr.  DBlhfeldf  B  Director , 

Deputy  for  Development  Plans,  ESD-XR,  (617)-271-3832) . 
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B.  Objectives 

Virtually  all  research  and  development  in  the  area  of  adaptive 
arrays  has  dealt  with  signal-to-noise  ratio  improvement  by  minimization 
of  directional  interference  and  ambient  noise  using  adaptive  signal 
processing.  Interference  of  this  nature  is  not  correlated  with  the 
desired  signal  and,  as  such,  can  be  called  "signal-uncorrelated" 
interference.  Another  equally  important  mechanism  of  signal  corruption 
is  the  linear  distortion  of  desired  signals  caused  by  reflections, 
scattering,  and  multipath.  In  wideband  digital  communication  systems, 
for  example,  this  introduces  intersymbol  interference,  detection 
errors  and,  consequently,  significant  degradation  in  performance.  Each 
array  element  or  each  array  beam  is  connected  to  an  adjustable  filter 
or  equalizer  (a  transversal  filter  or  tapped  delay  line)  and  each 
individual  filter  must  be  adjusted  in  a  systematic  way  to  minimize 
signal-correlated  interference.  The  research  effort  reported  here  is 
concerned  with  the  specific  problem  of  equalizer  adjustment  algorithms 
and  their  performance.  The  equalization  problem  is  formulated  in  the 
time  domain  as  an  integral  equation  of  convolution  which  is  solved 
in  the  minimax  sense  using  the  Remez  algorithm.  The  discrete-time 
case  is  also  considered  and  solved  via  the  Ascent  algorithm. 
Least-square-error  and  envelope-constrained  algorithms  also  are 
included. 

C.  Status  of  Research 

The  fundamental  research  problem  is  equalizer  design  in  the  til 
domain.  The  minimax  (Chebyshev)  error  criterion  was  chosen  so  that 
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equalization  errors  would  be  both  minimized  and  bounded  in  the  time- 
domain.  Since  frequency-domain  techniques  are  not  used,  hardware  and 
software  complexity  and  cost  may  be  reduced  and  filter  adjustment 
speeds  increased.  The  studies  reported  here  indicate  that  minimax 
time-domain  design  of  an  equalizer  is  both  reasonab’-i  and  practical. 
Experimental  results  from  computer-aided  designs  of  minimax  equalizers 
illustrate  several  aspects  of  the  performance  capabilities  and  compu¬ 
tational  sensitivities  of  the  Remez  algorithm.  Of  particular  interest 
are  the  trades  between  filter  length,  total  computation  time,  and 
accuracy.  Several  representative  examples  are  given.  Equations 
describing  all  filter  adjustment  algorithms  are  included.  Computer 
codes  for  the  Remez  algorithm  and  Ascent  algorithm  are  listed. 


D.  Technical  Publications  Related  to  Research 

1.  D.  Preis  and  C.  Bunks,  "Minimax  Time-Domain  Deconvolution 
for  Transversal  Filter  Equalizers,"  ICASSP-80  Proceedings, 
pp.  943-946,  IEEE  International  Conference  on  Acoustics, 
Speech  and  Signal  Processing,  Denver,  Colorado,  April  1980 

2.  D.  Preis,  and  C.  Bunks,  "Chebyshev  Time-Domain  Deconvolution 
for  Transversal  Filter  Equalizers,"  to  appear  in  Electronics 
Letters,  December  1980 

3.  C.  Bunks,  "Minimax  Equalization  for  Transversal  Filters," 
Master  of  Science  Thesis,  Tufts  University,  101  pages, 

June  1980 

4.  D.  Preis  and  C.  Bunks,  "Miniraax  Time-Domain  Deconvolution 
for  Transversal-Filter  Equalizers,  presented  at  the  IEEE 
International  Conference  on  Acoustics,  Speech  and  Signal 
Processing,  Paper  DSP  11.1,  Fairmont  Hotel,  Denver, 

Colorado,  April  1980 

5.  D.  Preis,  "Minimax  Equalizers:  Performance  and  Adjustment 

j  Algorithms,"  presented  at  the  1980  L’Aquila  Workshop  on 

‘  Digital  Signal  Processing,  sponsored  by  IEEE  Acoustics, 

Speech  and  Signal  Processing  Society,  L’Aquila,  Italy, 
September  1980 
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6.  D.  Preis  and  C.  Bunks,  "Computational  and  Performance  Aspects 
of  Minimax  Equalizers,"  accepted  for  presentation  at  IEEE 
International  Conference  on  Acoustics,  Speech  and  Signal 
Processing,  Sheraton-Atlanta  Hotel,  Atlanta,  Georgia,  March  1981 

7.  D.  Preis  and  C.  Bunks,  "Three  Algorithms  for  the  Design  of 
Transversal-Filter  Equalizers,"  submitted  for  presentation  at 
IEEE  International  Symposium  on  Circuits  and  Systems,  Chicago, 
Illinois,  April  1981 

E.  Professional  Personnel 

Biographical  Sketch  of  Principal  Investigator.  Douglas  Preis 
received  the  B.S.E.E.  degree  in  1964  and  M.S.E.E.  degree  in  1966,  both 
from  the  University  of  Santa  Clara,  CA.  He  was  a  NASA  Fellow  at  Utah 
State  University,  Logan,  where  he  received  the  Ph.D.  degree  in  1969. 

He  joined  the  research  faculty  of  Harvard  University,  Cambridge,  MA, 
in  1970  where  he  was  Post-Doctoral  Research  Fellow  in  the  Division  of 
Engineering  and  Applied  Physics  at  the  Gordon  McKay  laboratory  until 
1978.  In  1978,  he  joined  the  faculty  of  Tufts  University,  Medford,  MA, 
as  Assistant  Professor  of  Electrical  Engineering.  Dr.  Preis  has 
published  or  presented  about  thirty-five  original  technical  papers  and 
reports  in  the  electrical  engineering  field.  He  received  the  Audio 
Engineering  Society  Publications  Award  for  his  paper,  "Linear 
Distortion,"  which  was  published  in  the  Journal  of  Audio  Engineering 
Society  in  1976.  Dr.  Preis  is  a  member  of  Sigma  Xi,  Tau  Beta  Pi.  He 
is  a  member  of  the  Review  Board  of  the  Journal  of  the  Audio 
Engineering  Society  and  an  Associate  Editor  of  the  IEEE  Transactions 
on  Acoustics ,  Speech,  and  Signal  Processing. 

Biographical  Sketch  of  Research  Assistant.  Carey  Bunks 
simultaneously  received  the  B.S.  degree  summa  cum  laude  and  M.S. 
degree  from  Tufts  University  both  in  electrical  engineering  in  June  1980. 
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He  received  the  Amos  Emerson  Dolbear  Scholarship  and  the  Victor  Prather 
Prize  for  demonstrated  excellence  in  scientific  research  while  at  Tufts. 
He  presented  a  technical  paper  at  ICASSP-80  in  Denver  and  is  co-author 
of  four  other  technical  papers.  He  is  currently  pursuing  the  Ph.D. 
degree  in  electrical  engineering  at  MIT.  Mr.  Bunks  is  a  member  of  the 
Tau  Beta  Pi  and  Eta  Kappa  Nu. 
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MINIMAX  TIME-DOMAIN  DECONVOLUTION 
FOR  TRANSVERSAL  FILTER  EQUALIZERS 


C.  Bunks  and  D.  Preis 


Tufts  University 

Department  of  Electrical  Engineering 
Medford,  Massachusetts  02155 


ABSTRACT 


TIME-DOMAIN  DECONVOLUTION 


The  subject  of  this  paper  is  the  design  of  an 
optimum  transversal  filter  for  time-domain 
equalization  of  a  linear,  time-invariant  system. 
Given  advance  specification,  in  the  time  domain, 
of  the  unequalized  system  response  and  the  desired 
equalized  response,  a  numerical  deconvolution 
Is  performed  to  determine  the  set  of  filter  cap 
weights  which  optimally  satisfies  the  continuous¬ 
time,  convolutional  Integral  equation  in  the 
minimax  or  Chebyshev  sense.  This  solution  Is 
computed  using  the  second  algorithm  of  Remez. 
Numerical  examples  are  provided  which  illustrate 
the  efficacy  of  minimax  deconvolution.  A  computer 
program  flow  chart  is  included. 


INTRODUCTION 

Deconvolution  arises  In  practical  engi¬ 
neering  problems  as  diverse  as  time-domain 
equalization.  Instrumentation  and  measurement, 
system  Identification,  ultrasonic  diagnostics, 
geophysical  exploration,  and  imaging.  Convolution 
Is  itself  an  Integral  operation  whereas  deconvolu¬ 
tion,  or  the  Inverse  of  convolution,  requires  that 
an  Integral  equation  be  solved  either  directly  or 
Indirectly.  In  most  linear  signal  processing 
applications,  deconvolution  Is  accomplished 
numerically  and  approximately  using  either  itera¬ 
tive  techniques  flj  or  discrete  Fourier  transform 
methods  [2].  When  exact  solutions  do  not  exist, 
the  usual  procedure  la  to  3eek  an  approximate 
solution  having  minimum  square  error  [3].  Large 
errors  as  well  as  non-physical  results  sometimes 
associated  with  this  kind  of  solution  can  be 
suppressed  by  Imposing  suitable  constraints  when 
deconvolving  [l],  [4],  An  alternative  approach 
is  a  mlnlmax  or  Chebyshev  approximate  solution  [5] 
wherein  maximum  individual  errors  are  minimized 
rather  than  the  cumulative  square  error.  In  this 
report,  time-domain  equalization  of  a  continuous- 
time  system  with  an  N  tap  transversal  filter  is 
considered  as  a  deconvolution  problem.  The 
convolutional  integral  equation  Is  solved  In  Che 
Chebyshev  or  mlnlmax  sense  using  the  second 
algorithm  of  Remez  [5].  Maximum  deconvolution 
errors  are  controlled  by  the  number  of  taps  and 
minimized  uniformly. 


The  overall  Impulse  response  of  a  linear 
system  In  cascade  with  an  equalizing  filter  is 
determined  by  convolution  in  the  time  domain. 
Assuming  that  the  system  to  be  equalized  is 
causal  with  Impulse  response  h(t),  and  that  there 
is  a  particular  desired  response  g(t),  then  the 
filter’s  impulse  response  fit)  must  satisfy, 

f  f(x)h(t-x)dx  “  g(t).  (1) 

o  ' 

Since  the  unknown  filter  response  appears  under 
the  Integral  sign,  (1)  Is  an  integral  equation 
for  fit).  The  solution  of  (1)  requires  that  an 
Inverse  convolution  or  deconvolution  be  performed 
to  find  £(t)  when  h(t)  and  g(t)  are  specified. 


Given  a  transversal  filter  with  N  tap  weights, 
assume  that  the  filter's  finite-duration  impulse 
response  can  be  represented  by. 


fit) 


fn«(t-nT), 


(2) 


where  fn  represents  the  value  of  the  n-th  tap 
weight,  6  is  the  unit  impulse  function,  and  T  is 
the  time  delay  between  the  equally-spaced  taps. 

(An  alternative  representation  for  f(t)  is  a  set 
of  N  contiguous  rectangular  pulses  with  ampli¬ 
tudes  fn  and  width  T.  The  chosen  representation 
simplifies  the  following  development.) 

Substituting  (2)  into  integral  equation  (1)  yields 
the  following  continuous- time  representation  for 
8(t), 

N 

f  h(t-nT)  -  g(t).  (3) 

n 

This  result  implies  Chat  g(t)  can  be  represented 
exactly  by  a  superposition  of  N  weighted  and 
time-shifted  versions  of  h(t).  Depending  on 
both  hit)  and  the  choice  of  g(t).  It  may  not  be 
possible  to  satisfy  (3)  exactly.  In  such  cases, 
the  solution  to  this  problem  requires  finding  a 
set  of  N  tap  weights  fn  which  satisfies  (3) 
approximately.  For  any  specific  Instant  in  time  ta 
within  the  domain  of  definition  of  g(t),  13)  will 
yield  the  following  linear  equation  with  the  N  tap 
weights  as  unknowns, 


n  «  1 
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(4) 


z 


f  h(t  -nT) 
n  m 


g(t ). 
m 


If  M  distincc  Instants  in  time  are  chosen,  then  a 
system  of  M  linear  equations  In  N  unknowns  Is 
generated  which  takes  the  form, 

N 

f  h(t-nT)  -  g(t  ), 
n  in  in 

n  *  1 

m  -  1 ,  2,  ...  M.  (5) 

Since  It  may  not  be  possible  to  solve  (3)  exactly, 
a  systematic  method  will  be  developed  to  find  an 
approximate  solution.  This  is  accomplished  first 
by  considering  solutions  to  (3). 

Equation  (5)  suggests  three  possibilities: 

M  <  N,  M  «  N,  and  M  >  N.  For  M  <  N,  there  are 
more  unknowns  than  equations  so  chat  the  system 
of  linear  equations  is  underspecified  and,  con¬ 
sequently,  an  infinite  number  of  solutions  exists. 
For  M  «  N,  the  number  of  equations  equals  the 
number  of  unknowns  and  a  unique  solution  exists  if 
these  equations  are  linearly  independent.  This 
approximate  solution,  known  as  collocation  [6]  or 
point  matctiing,  only  satisfies  (3)  at  each  of  the 
chosen  N  points.  A  different  selection  of  N  points 
generally  will  yield  a  different  solution.  It  is 
clear  in  this  situation  that  the  tap  weights  fn 
depend  upon  the  choice  of  the  N  values  of  t_. 
Furthermore,  even  though  the  solution  to  (3) 
matches  exactly  at  the  selected  N  points  tn,  the 
approximation  error  may  be  considerable  between 
these  points.  For  M  >  N,  the  system  of  linear 
equations  is  overspecified  because  there  are  more 
equations  than  unknowns.  In  such  an  overspecified 
system  no  exact  solution  exists  which  simulta¬ 
neously  satisfies  each  of  the  M  equations.  Here 
the  only  recourse  is  to  find  the  "best"  approximate 
solution  to  (5).  A  popular  method  for  solving  (5) 
approximately,  which  also  ensures  a  minimum  total 
approximation  error,  is  the  method  of  least  squares 
[6]  or  minimum-square-error  technique.  In  this 
method  the  sum  of  the  square  of  each  of  the 
individual  errors  between  the  approximation  and 
the  desired  solution  of  the  M  equations  is 
minimized.  However,  in  spite  of  the  fact  that 
the  cumulative  error  for  the  total  system  of 
M  equations  is  minimized,  the  error  for  any 
Individual  equation  is  unconstrained  and  may  be 
very  /.urge  Indeed.  Therefore,  it  is  desirable  to 
seek  an  approximate  solution  which  simultaneously 
and  uniformly  minimizes  the  maximum  error  from  each 
of  the  M  individual  equations.  The  Chebyshev 
solution  to  an  overspecified  system  of  linear 
equations  satisfies  these  requirements  (3] . 

In  order  to  formulate  the  Chebyshev  solution 
to  (5),  it  is  necessary  to  define  an  error  equation 
associated  with  each  of  the  M  linear  equations 
having  the  general  form. 


r<t  > 

in 


N 

g(t  )  -  V'  fh(t-nT), 
m  n  ® 

n  *  1 

1,  2,  ...  M. 


(6) 


Next,  the  objective  is  to  minimize  the  maximum 
value  of  the  magnitude  of  r(tm)  for  a  ■  1,  2,  ...  M, 
by  appropriate  adjustment  of  the  N  tap  weights  fn. 
The  domain  of  definition  of  a  continuous  function 
such  as  g(t)  is  an  interval  containing  an  infinite 
number  of  points.  This  fact  implies  that  (3)  can 
be  viewed  as  an  overspecified  system  of  M  linear 
equations  in  N  unknowns  in  the  limit  as  M  -*  •  . 

Thus,  there  is  a  corresponding  infinite  number  of 
error  equations,  as  specified  in  (6),  which  define 
an  error  function  r(t)  of  the  form. 


r  (t) 


g(t)  - 


N 


f  h(t-nT) 
n 


(7) 


From  a  different  viewpoint,  the  location  of  each 
tm  (for  m  «  1,  2,  ...  M  and  M > N)  remains 
unspecified  yet  it  is  required  that  the  magnitude 
of  each  corresponding  individual  error  r(tn)  be 
minimized.  So,  it  is  necessary  to  examine  r(tm) 
for  all  values  of  t„  within  the  domain  of  defini- 
tion  of  g(t).  Specifically,  the  maximum  magnitude 
of  r(t)  must  be  minimized.  This  can  be  accom¬ 
plished  by  Invoking  the  characterization  theorem 
which  predicts  that  the  maximum  value  of  each  of 
the  M  >  N  error  equations  is  minimized  at  the 
intersection  of  a  particular  subset  of  N  +  1  of 
these  M  equations  [5].  Therefore,  finding  the 
Chebyshev  minimum  of  the  set  of  error  equations 
reduces  to  determining  the  appropriate  subset  of 
N  +  1  error  equations,  or,  equivalently,  finding 
the  best  3et  of  N  +  1  points  in  time  (viz.,  each 
tm  for  m  “  1,  2,  ...  N+l).  A  major  difficulty 
here  is  that  for  M  infinite,  an  infinite  number 
of  subsets  each  consisting  of  N  +  1  points  must  be 
examined.  A  systematic  search  for  the  subset  of 
N+l  points,  from  the  infinite  set  of  points  on 
the  domain  of  definition  of  g(t),  is  known  as  the 
second  algorithm  of  Remez  [5].  It  is  a  powerful 
iterative  technique  which  converges  to  the  set  of 
N+l  points  for  which  error  function  r(t)  is 
minimized  in  the  Chebyshev  sense.  The  N  tap 
weights  thereby  obtained  yield  the  requisite  mini¬ 
max  solution  to  (3).  This  entire  procedure  is 
defined  as  minimax  deconvolution.  It  is  the  best 
way  to  minimize  deconvolution  errors  uniformly. 

COMPUTED  EXAMPLES 


Figure  1  illustrates  equalization  of  a 
Gaussian  to  approximate  a  sine  function  on  the 
time  Interval  [-3*,  3s]  first  using  a  6  tap 
(T  -  it)  then  an  8  tap  (T  *  3»/4)  transversal 
filter.  In  this  example,  h(C)  ■  0.337exp(-t^/27.6) 
and  g(t)  ■  (sin  t)/t.  Maximum  value  of  the  error 
magnitude  |r(t)|  is  0.21  for  6  taps  and  is  reduced 
(by  a  factor  of  ten)  to  0.021  for  8  taps. 
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Figure  2  shows  a  raised-cosine  pulse 
equalized  Co  approximate  a  narrow  triangular  pulse 
on  the-  time  Interval  [-*,*]  using  an  6  tap 
(T  •  it/ 4)  then  a  12  tap  (T  -  ir/6)  transversal 


filter.  For  this  second  example,  h(t)  «  0.5  x 
(1  +  cos  t)  -  cos2(t/2);  g(t)  -  2  -  4|t|/x  for 
|t|  <  it / 2  and  g(t)  =  0  for  |t|  >  ir/2.  Maximum 
value  of  the  error  magnitude  |r(t)|  is  0.263  for 
8  taps  and  is  reduced  to  0.159  for  12  taps.  On  a 
percentage  basis,  these  errors  are  13%  and  8%, 
respectively,  because  the  peak  value  of  g(t)  is  2. 

COMPUTER  PROGRAM 

The  computer  program  (see  Figure  3)  first 
selects  N  +  1  evenly-spaced  points  t  in  the 
solution  interval  [a,b].  Each  of  these  points 
determines  one  linear  equation  with  the  N  tap 
weights  as  unknowns.  The  minimax  solution  to  this 
overspecified  system  of  N  +  1  equations  in  N 
unknowns  is  computed.  Next,  the  resulting 
continuous-time  error  function  r(t)  is  evaluated. 
If  the  maxima  of  |r(t)|  do  not  occur  at  the  chosen 
values  of  t  then  the  continuous-time  minimax 
solution  has  not  been  found.  The  next  iteration 
begins  by  selecting  new  points  tn  in  the  interval 
[a,b]  which  correspond  to  the  N  +  1  largest  values 
of  |r(t)[.  The  minimax  solution  to  this  new 
overspecified  system  is  computed  and  its  error 
function  r(t)  evaluated.  The  process  is  repeated 
until  | r (t ) |  is  bounded  by  |r(tm)|.  When  this 
occurs,  the  appropriate  subset  of  S  +  1  points  t 
is  known  and  the  minimax  solution  to  its  cor¬ 
responding  overspecified  system  of  N  +  1  linear 
equations  in  N  unknowns  yields  that  set  of  N  tap 
weights  which  satisfies  the  convolutional  integral 
equation  in  the  minimax  sense.  In  summary,  the 
Remez  algorithm  is  a  systematic,  iterative 
procedure  for  finding  a  unique  set  of  8  +  1  points 
tm  in  the  interval  exchanging  one  set  of  N  +  1 
points  tm  for  another  until  each  point  corresponds 
to  an  extreme  value  of  the  error  function  r(t). 
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-  Abstract  - 

Minimax  time-domain  deconvolution  to  determine  the  N  tap  weights 
of  a  transversal  filter  equalizer  (FIR)  is  accomplished  iteratively 
using  the  Remez  algorithm.  This  minimax  design  can  be  converted  to  a 
recursive  (IIR)  filter  possibly  having  greater  error  but  only  half 
the  number  of  delay  elements  T.  Results  of  computations  for  three 
different  examples  illustrate  the  dependence  of  equalization  error  e 
on  the  number  of  taps  N  and  width  of  solution  interval  W. 

Computation  time  and  convergence  are  related  to  N.  Representative 
FIR  and  IIR  designs  are  compared. 
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INTRODUCTION 


The  subject  of  this  contribution  is  the  automatic  design  of  an 
optimum  transversal  filter  for  time-domain  equalization  of  a  linear 
system.  The  equalization  problem  is  formulated  in  the  time  domain  as 
an  integral  equation  of  convolution.  Frequency-domain  techniques  are 
not  used.  Given  advance  specification  (or  an  estimate),  in  the  time 
domain,  of  the  unequalized  system  response  and  the  desired  equalized 
response,  a  numerical  deconvolution  is  performed  to  determine  the  set 
of  filter  tap  weights  which  optimally  satisfies  the  convolutional 
integral  equation  in  the  minimax  or  Chebyshev  sense.  The  transversal- 
filter  equalizer  is  assumed  to  be  a  tapped  delay  line  having  either  a 
finite-duration  impulse  response  or  an  infinite-duration  impulse 
response  while  the  linear  system  to  be  equalized  is  continuous- time. 

The  filter  tap  weights  are  computed  interatively  using  the  second 
algorithm  of  Remez.  A  major  advantage  of  minimax  time-domain  deconvolu¬ 
tion  for  equalizer  design  is  that  the  error  magnitude  is  both  minimized 
and  bounded.  The  question  of  solution  feasibility,  which  arises  in 
envelope-constrained  deconvolution,  is  avoided.  Also,  the  problem  of 
large,  unconstrained  individual  errors,  which  often  occurs  with 
least-square-error  approximate  solutions  or  discrete-Fourier-transform 
solutions,  is  eliminated. 

Specific  aspects  of  the  equalization  problem  also  considered  here 
Include  an  inter comparison  of  the  performance  of  both  recursive  and 
non-resursive  rainimax  transversal-filter  equalizers  in  terms  of 
equalization  accuracy  and  filter  length. 
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Prospective  applications  of  this  research  include:  high-data- 
rate  communication  through  time-varying  dispersive  channels,  reduction 
of  intersymbol  Interference,  system  identification  and  modeling,  minimax 
inverse  filtering,  multichannel  filtering,  and  sensor  array  equalization. 
Potential  advantages  of  minimax  time-domain  equalization  are:  a  reduc¬ 
tion  of  hardware  and  software  complexity  and  costs  normally  associated 
with  time-domain  to  frequency-domain  and  frequency-domain  to  time-domain 
transformations,  increased  processing  speed,  uniform  control  and  mini¬ 
mization  of  deconvolution  errors  in  the  time  domain,  and  the  possibility 
of  reduced  filter  complexity. 

Computed  results  for  about  50  FIR  equalizer  designs  are  summarized 
in  graphical  forms  which  relate  accuracy,  convergence,  and  computation 
time  to  filter  length.  Also  illustrated  are  sensitivities  of  equaliza¬ 
tion  error  to  solution- interval  width.  These  minimax  designs  were 
computed  using  the  Remez  algorithm  [1]. 

EXAMPLES  AND  RESULTS 

Figure  l  summarizes  three  different  time-domain  equalizer  design 
problems.  In  the  left  column  of  this  figure  are  the  key  symbols  which 
are  used  in  subsequent  figures  to  identify  each  of  these  three 
examples.  The  first  case  is  equalization  of  a  Gaussian  to  a  Sine 
function,  the  second,  a  raised  cosine  to  a  triangular  function  and 
the  third,  a  raised  cosine  to  a  rectangular  function.  The  filter  or 
equalizer  impulse  response  is  f(t),  the  system  response  h(t),  and 
the  overall  response  g(t).  The  minimax  equalizer  is  a  tapped  delay 
line  or  transversal  filter  having  finite-duration  impulse  response  (FIR). 


In  some  examples  to  follow  infinite  duration  impulse  response  (IRR) 
filters  also  will  be  considered.  Note  that  h(t)  and  g(t)  are  continuous¬ 
time  functions  while  the  filter  f(t)  is  discrete-time  device. 

Two  equalizer  canonical  forms  are  shown  in  Figure  2.  The  letter 
T  indicates  a  time  delay  of  T  seconds.  Lower  case  letters  f,  a,  and  b 
designate  filter  coefficients  or  tap  weights. 

Figure  3  illustrates  the  dependence  of  equalization  accuracy  on  the 
number  of  taps  N  for  the  FIR  designs.  The  number  of  taps  equals  the 
number  of  delay  elements  minus  one.  Since  the  filter  output  is  a  weighted 
sum  of  time  delayed  versions  of  the  input,  the  symmetry  of  the  taps 
about  t  3  0  is  important.  For  N  odd  equalization  accuracy  is  better  than 
N  even  because  the  middle  tap  is  centrally  located  at  t  »  0.  Considering 
N  even  or  N  odd  separately  it  is  seen  that,  for  fixed  solution  interval 
width  W,  the  error  decreases  the  equalization  accuracy  improves  as  the 
number  of  taps  N  is  increased.  The  error  is  computed  as  the  ratio  of 
maximum  deviation  of  equalized  response  from  g(t)  to  maximum  value  of 
g(t)  and  this  is  expressed  in  percent.  Note  that  for  the  third  example 
(square  symbol)  the  error  can  not  be  less  than  50%  because  (at  best) 
the  equalized  response  will  pass  through  the  vertical  sides  of  the 
rectangle  at  their  midpoints. 

In  Figure  U  minimax  equalization  error  is  plotted  for  the  three 
test  cases  where  the  number  of  taps  is  fixed  at  N  *  8  and  the  width 
of  the  solution  interval  W  is  varied.  It  seems  intuitive  that  the 
error  should  increase  as  the  solution  interval  is  made  wider, 
however,  this  is  not  necessarily  true  as  can  be  observed  for  the 
first  example  (round  dots).  Also,  it  is  interesting  to  note  that  if 
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the  solution  interval  is  chosen  to  be  wider  than  the  sum  of  the 
durations  of  h(t)  and  f(t)  then  no  change  in  error  can  be  expected 
since  the  convolution  is  zero  for  |t|  large. 

Convergence  of  the  Remez  algorithm  depends  on  several  factors. 

Figure  5  illustrates  convergence  as  a  function  of  the  number  of  tap 
weights  N.  It  is  seen  that  filter  symmetry  (N  odd  or  even)  is  influential. 
On  the  average,  the  number  of  iterations  required  for  convergence  is 
approximately  equal  to  the  number  of  taps  N. 

Figure  6  shows  total  computation  time  in  seconds  as  a  function  of 
the  number  of  tap  weights  N.  With  one  exception,  computation  time 
increases  with  N.  Those  filters  offering  higher  accuracy  (N  odd) 
require  proportionally  more  computation  time. 

Figure  7  combines  the  data  of  Figures  5  and  6  to  illustrate 
the  average  time  per  iteration  as  a  function  of  N.  Note  that  the 
increased  time  for  the  first  example  (round  dots)  is  due  to  the  compara¬ 
tively  more  complicated  form  of  the  mathematical  expression  for  g(.t) 
and  h(t) . 

Figure  8  shows  the  trade-off  between  filter  accuracy  and  invested 
computation  time.  Equalization  error  can  be  reduced  significantly 
by  lengthening  the  filter.  For  short  filters,  computation  time  and 
accuracy  are  related  linearly. 

Figure  9  is  a  table  of  two  sets  of  computed  tap  weights  for 
nonrecursive  (FIR)  and  recursive  (HR)  designs  when  N  ■  7,  T  =  ir/7, 
for  W  *  6ir,  and  for  W  ■  12ir.  Tap  weights  for  the  IIR  case  were 
computed  from  the  FIR  values  using  an  algorithm  [2]  which  forces  the 
initial  N  values  of  the  impulse  response  of  the  IIR  filter  to  equal 
;he  corresponding  values  of  the  FIR  filter. 
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Figure  10  compares  FIR  and  HR  equalized  responses  with  g(jt). 

The  maximum  FIR  error  in  the  solution  interval  W  ■  6ir  is  1.52%  for 
N  =  7. 

Figure  11  illustrates  the  large  errors  (side  lobes)  outside  the 
solution  interval  W  =  6ir  for  the  designs  of  Figure  10.  These  sidelobes 
can  be  reduced  by  increasing  W  as  shown  in  Figure  12  where  the  maximum 
error  is  13.3%. 
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Figure  2.  Equalizer  canonical  forms 
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Figure  3.  Dependence  of  minimax  error  on  number  of  taps 
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Figure  4.  Dependence  of  minimax  error  on  width  of  solution  interval 


CONVERGENCE 


Number  of  iteratons  for  convergence  vs.  number  of  taps 


Figure  6.  Computation  time  vs.  number  of  taps 


AVERA 


Average  time  per  iteration  vs.  number  of 
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Figure  9.  Tap  weights  for  FIR  and  HR  equalizer  using  in  figures 
to  follow 
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6tt  solution  interval  W. 
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Figure  11.  Same  as  previous  figure  but  including  response  outside 
solution  interval  W. 
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12tt  solution  interval 


CHEBYSHEV  TIME-DOMAIN 


DECONVOLUTION  FOR  TRANSVERSAL- 
FILTER  EQUALIZERS 


D.  Preis  and  C.  Bunks* 


-  Abstract  - 

Design  equations  are  presented  for  computing  the  N  tap  weights 
of  a  transversal  filter  equalizer.  The  associated  time-domain 
deconvolution  is  performed  subject  to  the  Chebyshev  (minimax) 
error  criterion.  Individual  or  pointwise  deconvolution  errors 
are  controlled  and  minimized  uniformly.  This  method  offers  an 
alternative  to  least-square-error  or  discrete-Fourier-transform 
designs  wherein  only  the  cumulative  deconvolution  error  is 
minimized.  Computationally,  the  Ascent  algorithm  is  used  which 
requires  inversion  of  two  N  +  1  by  N  +  1  matrices  followed  by  a 
finite  number  of  elementary  algebraic  exchange  operations.  The 
equations  needed  to  implement  this  algorithm  are  presented  in  a 
form  suitable  for  digital  computation. 

To  appear  in  Electronics  Letters 

*The  authors  are  with  the  Department  of  Electrical  Engineering, 
Tufts  University,  Medford,  Massachusetts  02155  USA. 
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When  a  transversal  filter  (tapped  delay  line)  is  used  as  an 

equalizer  for  a  discrete- time  system,  it  is  necessary  to  compute  the 

N  tap  weights  of  the  filter  by  deconvolution.  Because  the  operation 

of  discrete-time  convolution  generates  more  equations  than  unknowns, 

deconvolution  generally  will  be  inexact  and,  consequently,  equalization 

will  be  approximate.  Any  approximate  solution  to  this  equalizer  design 

problem  will  depend  on  the  error  criterion  and  constraints  imposed  when 

1  2 

deconvolving.  In  two  previous  communications  *  design  equations  were 
given  for  least-square-error  and  envelope-constrained  transversal- 
filter  equalizers.  The  present  design  procedure  offers  an  alternative 
wherein  the  maximum  deconvolution  errors  are  minimized  uniformly. 

This  is  known  as  the  Chebyshev  (or  minimax)  approximate  solution. 

Consider  the  discrete-time  description  of  equalization  given  by 
the  convolution  summation, 

g(m)  »  l  h  (m-n+1)  f(n),  (1) 

n=l 

where  g  is  a  specified  M  element  sequence  representing  the  desired 
equalized  Impulse  response,  f  is  an  N  element  sequence  composed  of 
the  unknown  transversal-filter  tap  weights,  and  h  is  a  known  M-.M+l 
element  sequence  which  represents  the  impulse  response  of  the  linear, 
time- invariant  system  requiring  equalization.  It  is  assumed,  for 
convenience,  that  the  time  interval  between  adjacent  sequence  elements 
is  unity.  In  all  practical  cases,  the  sequence  h  has  more  than  one 
element  so  that  M>N  and  eqn.  1  generates  an  overspecified  system 
of  M  linear  equations  in  N  unknowns.  Since  an  overspecified  system 
such  as  this  may  not  have  an  exact  solution,  it  is  appropriate  to 
consider  approximate  solutions.  A  useful  measure  of  the  approximation 
error  (or  equalization  error)  corresponding  to  a  specific  set  of  N 
tap  weights  f  is  the  M  element  error  sequence  r  defined  by 
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r(m)  -  g(m)  -  £  h(m-n+l)  f(n)  (2) 

n=l 

Comparison  of  eqn.  2  with  eqn.  1  will  indicate  that  the  closer  each 
numerical  value  r(m)  is  to  zero,  the  better  the  approximate  solution 
becomes.  If,  for  example,  r(m)  =  0  for  m=  1,  2,  3,  ...  M  then  the 
N  tap  weights  satisfy  eqn.  1  exactly. 

The  intent  here  is  to  outline  a  procedure  for  computing  those 
N  tap  weights  f  which  both  minimize  and  bound  all  the  elements  of 
the  error  sequence  r  ,  that  is, 

|r(m)|<e  ,  (3) 

where  e  is  the  smallest  possible  non-negative  number.  This  solution 
is  known  as  the  Chebyshev  (or  minimum-maximum-error)  approximate 
solution.  There  are  N+l  unknowns  for  a  given  equalizer  design,  namely, 
the  N  tap  weights  and  e  .  Intuitively,  these  unknown  quantities  can  be 
found  by  solving  an  appropriate  set  of  N+l  linear  equations  in  N+l 
unknowns.  It  turns  out  that  this  set  is,  in  fact,  a  subset  of  the 
M  equations  specified  by  eqn.  1.  Assuming  e>0  the  Chebyshev  solution 
to  these  N+l  equations  is  characterized  by  the  following  two 


properties.  First, 


r(mi)  =  a(mi)E 


for  N+l  of  the  M  elements  of  the  error  seauence  given  by  eqn.  ? 

where  aOn^)  equals  +1  or  -1  and  the  subscript  i  =  1,  2,  3,  ...  N+l 

denotes  N+l  specific  integer  values  of  m  from  the  total  set  of  M 

values.  Second,  the  remaining  M  -  (N+l)  elements  of  r  will  satisfy 

eqn.  3.  One  systematic  search  for  the  appropriate  set  of  N+l  linear 

equations  and  the  subsequent  solution  for  the  N  tap  weights  f  and 

3 

maximum  error  e  is  known  as  the  Ascent  algorithm.  The  equations 


required  to  implement  this  algorithm  follow. 
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Initiate  the  search  by  selecting  any  subset  of  N+l  of  the  M 
equations  given  by  eqn.  2.  Using  eqn.  4  write  the  result  in  the 
following  matrix  form, 


j  g(m^)  o(m^)  h(m^)  h(m^-l)  .  .  .  h(m^-N+l) 

g(m2)  a(m2)  h<m2>  h(m2~l) 

g(m3>  a(m3)  h(m3)  hCn^-l) 


a(m2)  h(m2)  h(m2~l) 


.  h(m2-N+l) 
.  h(m3-N+l) 


[8(n,N+i)J  h(nW  h(mN+r1)-  •  •  h(raN+rN+1)  J  [f(N)J 

where  the  subscripts  associated  with  the  index  m  denote  those  specific 
IH-1  integer  values  of  m  selected  from  the  M  available  values.  Each 
of  the  N+l  values  of  a  appearing  in  eqn.  5  is  either  +1  or  -1  and  is 
specified  by  the  relation 

a(m1)  *  sgn  [0 (mi> ] ,  (6) 

where  the  numerical  value  of  each  of  the  N+l  elements  of  8  is  found 


by  inverting  the  following  matrix  equation, 


h(ra^)  h(m2)  h(m3> 

h(m;L-l)  h(m2-l)  h(m3-l) 


0(11^) 

0(m2) 


h^mN+l^  0(m2) 

hCn^+i"1)  0(m3) 


M  .  .  j 

h(mx-N+l)  h(ra2-N+l)h(ra3-N+l)  .  .  .  Mm^-N+l)!  |e(mN+1)j  ^  0  j 

The  existence  of  the  Chebyshev  solution  requires  that  zero  be  contained 
in  the  convex  hull  of  the  N+l  vectors  whose  N  components  are  the 
elements  of  the  N+l  columns  appearing  below  the  top  row  of  the  square 
matrix  in  eqn.  7.  This  condition  is  satisfied  if  a  linear  combination 
of  these  N+l  vectors  can  be  found  which  yields  the  null  (zero)  vector 
and  if  those  N+l  scalars  9,  used  to  form  the  linear  combination,  sum 
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to  unity.  Equation  7  is  a  mathematical  statement  of  these  two 
requirements.  The  algebraic  sign  of  each  of  the  N+l  scalars  9 
determines  whether  the  error  e,  associated  with  each  linear  equation 
in  eqn.  5,  is  positive  or  negative. 

Write  the  solution  to  eqn.  5  in  matrix  form  as. 


fc(0,0)  c(0,l)  .  .  .  c(0 , j )  .  .  .  c(0,N) 

r  8<ni) 

e 

lc(l.O)  c(l,l)  .  .  .  c(l,j)  .  .  .  c(l,N) 

• 

g(m2^ 

f(l) 

•  ♦  •  • 

•  •  •  • 

c(k,0)  c(k,l)  .  .  .  c(k,j)  .  .  .  c(k,N) 

•  •  •  • 

— 

• 

•  •  •  • 

•  •  •  • 

c(N ,0)  c(N,l)  .  .  .  c(N,j)  .  .  .  c(N,N) 

i^+i^ 

f(N) 

(8) 


The  matrix  multiplication  prescribed  by  eqn.  8  yields  the  initial  set 
of  numerical  values  of  the  N  tap  weights  f  and  error  e  .  Insert 
these  computed  tap  weights  in  eqn.  2  and  examine  each  element  of  the 
error  sequence  r  to  verify  that  |r(m)|£e  for  m  =  1,  2,  3,  . . .  M.  If 
this  is  true  then  eqn.  3  is  satisfied  and  the  Chebyshev  solution  has 
been  found.  No  further  computations  are  necessary. 

If  eqn.  3  is  not  satisfied  then  the  following  iterative  exchange 
procedure  is  required.  Select  that  specific  integer  m=p  for  which 
[r(m)|  is  maximum  and  let 

U  =  sgn  [r (p) ]  (9) 

Next,  calculate  the  N+l  elements  of  sequence  aO)  defined  by 

N 

a(j)  -  c(0, j)  +  l  h(p-k+l)  c  (k.j)  (10) 

k=l 

Now  select  the  integer  q  so  that 

b(q)  -  ya(q)/c(0,q)  (11) 

is  a  maximum.  Next  replace  the  N+l  element  column  c(k,  q)  in  eqn.  8 
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with  c'(k,q)  where 

c'(k,q)  =  c(k,q)/a(q),  (12) 

and  replace  all  other  elements  c(k,j)  in  eqn.  8  for  which  j  i*  q  with 
c'(k,j)  where 

c’(k,J)  =  c(k,j)  -  a(j)c’(k,q)  (13) 

Finally,  replace  g(q)  with  g(p)  in  eqn.  8.  Recompute  the  error  e, 
and  the  sequences  f  and  r  using  eqn.  8  and  eqn.  2.  As  before, 
verify  that  |r(m)!_^e  for  m  =  1,  2,  3,  ...  M.  If  eqn.  3  is  not 
satisfied  then  repeat  the  foregoing  exchange  procedure  by  selecting 
a  new  value  of  p  and  returning  to  eqn.  9. 

Note  that  although  H  -  (N+l)  exchange  operations  are  possible, 
the  Chebyshev  solution  usually  is  found  with  fewer  iterations  than 
this.  The  maximum  error  e  will  increase  (ascend)  with  each 
successive  iteration. 

In  certain  applications,  h  and  g  are  continuous-time  functions 

rather  than  finite  sequences  of  sampled  values  as  in  eqn.  1.  If  this 

is  the  case,  then  the  Chebyshev  solution  can  be  found  using  the  Remez 
3 

algorithm.  Mathematical  details  and  representative  examples  are 

A 

available  elsewhere. 
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-  Abstract  - 

This  report  contains  the  equations  needed  to  implement  three 
different  equalizer  design  algorithms.  The  equations  are  presented 
in  matrix  form  and  are  suitable  for  digital  computation.  The  first 
design,  which  is  based  on  the  Chebyshev  or  minimax  error  criterion, 
is  accomplished  using  an  iterative  procedure  known  as  the  Ascent 
algorithm.  The  second  design  is  a  direct,  least-square  minimization 
of  the  cumulative  equalization  error  by  fast  Toeplitz  matrix 
inversion.  The  third  design  technique  attempts  to  enforce  specific 
constraints  on  the  envelope  of  the  equalized  response.  If  an 
envelope-constrained  solution  is  feasible,  then  it  can  be  found 
iteratively  beginning  with  the  least-square-error  solution. 

Computationally,  the  least-square  minimization  is  the  fastest 
while  the  envelope-constrained  design  is  the  slowest.  The  minimax 
design  offers  the  best  error  control. 


Submitted  for  presentation  at  The  1981  International  Symposium  on 
Circuits  and  Systems,  Chicago,  Illinois,  April  27-29. 

Suggested  Session:  Computer-Aided  Design,  Active  Filters,  or 
Digital  Filters 


V-l 


-2- 


1.  Introduction.  When  a  transversal  filter  (tapped  delay  line)  is 
used  as  an  equalizer  for  a  discrete-time  system,  it  is  necessary  to 
compute  the  N  tap  weights  of  the  filter  by  deconvolution.  Because 
the  operation  of  discrete-time  convolution  generates  more  equations 
than  unknowns,  deconvolution  generally  will  be  inexact  and,  conse¬ 
quently,  equalization  will  be  approximate.  Approximate  solutions  to 
this  equalizer  design  problem  depend  on  both  the  error  criterion  and 
constraints  imposed  when  deconvolving.  This  report  presents  three 
different  deconvolution  algorithms  to  compute  the  N  tap  weights  of 

a  transversal-filter  equalizer:  (1)  minimax  design,  (2)  least- 
square-error  design,  and  (3)  envelope-constrained  design. 

2.  Minimax  Design.  The  basis  for  the  minimax  design  is  a  uniform 
minimization  of  maximum  deconvolution  errors.  Computationally,  the 
Ascent  algorithm  is  used  to  determine  the  N  tap  weights  of  the 
filter  and  this  requires  inversion  of  two  N+l  by  N+l  matrices 
followed  by  a  finite  number  of  elementary  exchange  operations. 

Consider  the  discrete-time  description  of  equalization  given 
by  the  convolutional  summation, 
m 

g(ra)  =  Y!  h  (m  -  n  +  1)  f  (n)  ,  (1) 

n  =  1 

m  =  1,  2,  3,  .  .  .  M 
n  *  1,  2, 3,  .  .  .  N 

where  g  is  a  specified  M  element  sequence  representing  the  desired 
equalized  impulse  response,  f  is  an  N  element  sequence  composed  of 
the  unknown  transversal- filter  tap  weights,  and  h  is  a  known  M  -  S  +  1 
element  sequence  which  represents  the  impulse  response  of  the  linear. 
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time- invariant  system  requiring  equalization.  It  is  assumed,  for 
convenience,  that  the  time  interval  between  adjacent  sequence  elements 
is  unity.  In  all  practical  cases,  the  sequence  h  has  more  than  one 
element  so  that  M  >  N  and  eqn.  1  generates  an  overspecified  system 
of  M  linear  equations  in  N  unknowns.  Since  an  overspecified  system 
such  as  this  may  not  have  an  exact  solution,  it  is  appropriate  to 
consider  approximate  solutions.  A  useful  measure  of  the  approximation 
error  (or  equalization  error)  corresponding  to  a  specific  set  of  N 
tap  weights  f  is  the  M  element  error  sequence  r  defined  by 

m 

r(m)  =  g(m)  -  ^  ^h  (m  -  n  +  1)  f(n),  (2) 

n  =  1  m=l,2,3,...M 

n  =  1,  2,  3,  .  .  .  N 

Comparison  of  eqn.  2  with  eqn.  1  will  indicate  that  the  closer  each 
numerical  value  r(m)  is  to  zero,  the  better  the  approximate  solution 
becomes.  If,  for  example,  r(m)  =  0  for  m  =  1,  2,  3,  .  .  .  M  then 
the  N  tap  weights  satisfy  eqn.  1  exactly. 

The  intent  here  is  to  outline  a  procedure  for  computing  those 
N  tap  weights  f  which  both  minimize  and  bound  all  the  elements  of  the 
error  sequence  r  ,  that  is , 

|r(m)  |_<  e,  (3) 

m  =  1,  2,  3,  .  .  .  M 

where  e  is  the  smallest  possible  non-negative  number.  This  solution 
is  known  as  the  Chebyshev  (or  miniraum-maximura-error)  approximate 
solution.  There  are  N  +  1  unknowns  for  a  given  equalizer  design, 
namely,  the  N  tap  weights  and  e  .  Intuitively,  these  unknown  quantities 
can  be  found  by  solving  an  appropriate  set  of  N  +  1  linear  equations  in 
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N  +  1  unknowns.  It  turns  out  that  this  set  is,  in  fact,  a  subset  of 
the  M  equations  specified  bv  eqn.  1.  Assuming  e  >0,  the  Chebyshev 
solution  to  these  N  +  1  equations  is  characterized  by  the  following 
two  properties.  First, 


rCnr)  =  (4) 

for  N  +  1  of  the  M  elements  of  the  error  sequence  given  by  eqn.  2, 

where  a(nu)  equals  +1  or  -1  and  the  subscript  i=l,  2,3,  .  .  .N  +  l 

denotes  N+l  specific  integer  values  of  m  from  the  total  set  of  M 

values.  Second,  the  remaining  M  -  (N+l)  elements  of  r  will  satisfy 

eqn.  3.  One  systematic  search  for  the  appropriate  set  of  N  +  1  linear 

equations  and  the  subsequent  solution  for  the  N  tap  weights  f  and 

1 

maximum  error  e  is  known  as  the  Ascent  algorithm.  The  equations 
required  to  implement  this  algorithm  follow. 

Initiate  the  search  by  selecting  any  subset  of  N  +  1  of  the  M 
equations  given  by  eqn.  2.  Using  eqn.  4  write  the  result  in  the 
following  matrix  form. 


1 

rH 

to 

1 _ 

a(m^)  h(m^)  h(m^-l)  .  .  .  h(m^-N+l) 

e 

g(m2) 

a(ra2>  h(m2)  h(m2-l)  .  .  .  h(m2~N+l) 

f(l) 

g(m3) 

= 

a  (m^)  h(m3)  hCm^-l)  .  .  .  h(m.j-N+l) 

•  •  •  • 

f(2) 

8{nW 

•  •  •  • 

•  •  ■  • 

°(Vi)h(nWi)  h(mN+r1)-  ■  •  h(mN+rN+1) 

f(N) 

where  the  subscripts  associated  with  the  index  ra  denote  those  specific 
’!  +  1  integer  values  of  m  selected  from  the  M  available  values. 
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Each  of  the  N  +  1  values  of  a  appearing  in  eqn.  5  is  either  +1  or  -1 
and  is  specified  by  the  relation 

=  sgn  [0(mi)3,  (6) 

i  =  1,  2,  3,  .  .  .  N  +  1 

where  the  numerical  value  of  each  of  the  N  +  1  elements  of  9  is  found 
by  inverting  the  following  matrix  equation. 


1 

l 

1 

...  1 

9  (mi) 

—  — 

1 

h(mj) 

h(m2) 

h(m3) 

.  .  .  hCm^) 

0(m2) 

0 

h(mj-l) 

• 

h(m2-l) 

h(m3~l) 

• 

•  •  • 

• 

0(m3) 

• 

= 

0 

• 

Mt^-n+i) 

h(m2-N+l) 

• 

h(m3-N+l) 

.  .  .  hCm^-N+l) 

• 

_9(,W 

0 

The  existence  of  the  Chebyshev  solution  requires  that  zero  be  contained 
in  the  convex  hull  of  the  N  +  1  vectors  whose  N  components  are  the 
elements  of  the  N  +  1  columns  appearing  below  the  top  row  of  the  square 
matrix  in  eqn.  7.  This  condition  is  satisfied  if  a  linear  combination 
of  these  N  +  1  vectors  can  be  found  which  yields  the  null  (zero)  vector 
and  if  those  N  +  1  scalars  0,  used  to  form  the  linear  combination,  sum 
to  unity.  Equation  7  is  a  mathematical  statement  of  these  two 
requirements.  The  algebraic  sign  of  each  of  the  N  +  1  scalars  9 
determines  whether  the  error  e,  associated  with  each  linear  equation 
in  eqn.  5,  is  positive  or  negative. 

Write  the  solution  to  eqn.  5  in  matrix  form  as. 
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c(0,0) 

c(0 ,1) 

.  .  .  c(0,j) 

.  .  .  c(0,N) 

gCn^) 

e 

c(l,0) 

c(l.l) 

.  .  .  c(l,j) 

• 

.  .  .  c(l,N) 

• 

g(m2) 

= 

f(l) 

c(k,0) 

c(k,l) 

• 

•  •  •  C  (k.  y  j  ) 

« 

• 

.  .  .  c(k,N) 

• 

• 

c(N  ,0) 

« 

1 

c(N,l) 

C 

ft 

.  .  .  c(N, j) 

.  .  .  c(N,N) 

g(mN+l) 

f(N) 

The  matrix  multiplication  prescribed  by  eqn.  8  yields  the  initial  set 
of  numerical  values  of  the  N  tap  weights  f  and  error  e  .  Insert  these 
computed  tap  weights  in  eqn.  2  and  examine  each  element  of  the  error 
sequence  r  to  verify  that  |r(m)  |_<  t  for  m  =  1,  2,  3,  .  .  .  M.  If 
this  is  true  then  eqn.  3  is  satisfied  and  the  Chebyshev  solution  has 
been  found.  No  further  computations  are  necessary. 

If  eqn.  3  is  not  satisfied  then  the  following  iterative  exchange 
procedure  is  required.  Select  that  specific  integer  ra  =  p  for  which 
[r(m)|  is  maximum  and  let 


y  =  sgn  [r (p) ] 


Next,  calculate  the  N  +  1  elements  of  sequence  a(j)  defined  by 


(9) 


a(j)  3  yc(0,j)  +  h(p  -  k  +  1)  c  (k,  j), 

j  -  0,  1,  2,  3,  . 


(10) 


k  =  1 


.  N. 


Now  select  the  integer  q  so  that 


1 


V-6 


-7- 


b(q) 


ua(q) 
c(0,  q) 


q  -  0,  1,  2,  3,  ...  N 


(ID 


is  a  maximum.  Next  replace  the  N  +  1  element  column  c(k,  q)  in  eqn.  8 
with  c'(k,  q)  where 


c' (k,q) 


-  c(k»  q) 

a(q)  ’ 


k  =  0,  1,  2,  3,  .  .  .  N 


(12) 


and  replace  all  other  elements  c(k,  j)  in  eqn.  8  for  which  j  4  q  with 
c'(k,  j)  where 


c'(k,  j)  =  c(k,  j)  -  a(j)c'(k,  q) ,  (13) 

k  -  0,  1,  2,  3,  ...  N 
j  =  0,  1,  2,  3,  .  .  .  N 

(j  *  q) 


Finally,  replace  g(q)  with  g(p)  in  eqn.  8.  Recompute  the  error  e, 
and  the  sequences  f  and  r  using  eqn.  8  and  eqn.  2.  As  before,  verify 
that  |r(m) (<_  e  for  m  =  1,  2,  3,  .  .  .  M.  If  eqn.  3  is  not  satisfied 
then  repeat  the  foregoing  exchange  procedure  by  selecting  a  new  value 
of  p  and  returning  to  eqn.  9. 

Note  that  although  M  -  (N  +  1)  exchange  operations  are  possible, 
the  Chebyshev  solution  usually  is  found  with  fewer  iterations  than 
this.  The  maximum  error  e  will  increase  (ascend)  with  each 
successive  iteration. 

In  certain  applications,  h  and  g  are  continuous-time  functions 
rather  than  finite  sequences  of  sampled  values  as  in  eqn.  1.  this 
is  the  case,  then  the  Chebyshev  solution  can  be  found  using  the  Remez 

algorithm. *  Mathematical  details  and  representative  examples  are 

2 

available  elsewhere. 
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3.  Least-Square-Error  Design.  The  least-square-error  design  is  based 
on  minimization  of  sum  of  the  squares  of  the  M  individual  deconvolution 
errors  r(m)  as  defined  in  eqn.  2.  Minimization  of  this  cumulative 
error  is  different  from  the  minimax  design  which  uniformly  and 
simultaneously  minimizes  each  individual  deconvolution  error  r(ra). 

Very  efficient  computation  of  the  least-square-error  solution  is 
possible  using  fast  Toeplitz  matrix  inversion. 

A  matrix  representation  of  the  convolutional  summation  given  in 
eqn.  1  is 


where  f  and  g  now  are  understood  to  be  column  vectors.  The  elements 
of  H  are  known  and  g  is  specified.  The  unknown  tap-weight  vector 
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is  f.  Thus,  eqn.  15  represents  an  overspecified  system  of  M  linear 

equations  in  N  unknowns,  that  is  M  >  N.  The  least-square-error 

-  3 

solution  f  to  eqn.  15  is 


/v 


f 


Sg 


(16) 


where  S  =  (  H'  H  )-1H  1 .  The  prime  symbol  denotes  transpose. 

Note  that  H'  H  is  an  N  x  N  symmetric  Toeplitz  matrix  whose  elements 

are  T. .  *  Ti.  .  i  =  T,  where 
ij  |i  “  J  I  k 


T 


k 


M-N+l-k 


h(n  +  k)  h(n) . 


(17) 


A  rapid  inversion  algorithm  is  available^  to  compute  ( H*  H) 

After  this  inverse  is  evaluated,  the  least-square-error  tap  weights 

/v 

f  are  found  using  eqn.  16.  The  resultant  mean-square  error  is 
[Hf  -  8]  [Hf  -  g]/M,  and  it  is  this  quantity  that  has  been  minimized. 
However,  the  individual  deconvolution  errors  r(m)  are  unconstrained. 
These  errors  can  be  evaluated  quantitatively  by  replacing  f(n)  with 

A  A 

f(n)  in  eqn.  2.  If  the  least-square-error  solution  f  yields 
unacceptably  large  individual  errors  r(m) ,  it  may  be  possible  to 
reduce  them  using  the  envelope-constrained  design  procedure  discussed 
in  the  following  section.  With  this  algorithm  somewhat  larger  mean- 

larger  individual  deconvolution  errors  r(m). 

A 

The  least-square-error  solution  f  also  can  be  found  using  discrete 
Fourier  transform  methods.  Computationally,  these  two  different 
methods  are  comparable.  One  requires  matrix  algebra  and  Toeplitz 
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inversion  while  the  other  relies  on  the  availability  of  FFT  routines 
and  requires  algebraic  operations  with  complex  numbers. 


4.  Envelope- Cons trained  Design.  The  envelope- cons trained  design 
algorithm  is  an  iterative  procedure  which  begins  with  the  least- 

/v 

square-error  design.  The  least-square-error  tap  weights  f  are  then 
systematically  re-adjusted  to  the  new  values  f  in  an  attempt  to  satisfy 
specific  constraints  imposed  on  the  envelope  of  the  equalized 
response  g(m). 

Assume  that  each  g(m)  is  bounded  above  by  a(m)  and  below  by  b(m). 
Next  require  that  the  envelope-constrained  tap  weights  f  satisfy 


Hf  - 


where 


b(m)  £  g(m)  <  a(m)  , 


m-1,  2,  3,  .  .  .  M. 


In  this  way,  the  equalized  response  g  is  constrained  to  lie  within 
the  tolerance  envelope  prescribed  by  sequences  a  and  b;  and  g  is  the 
envelope-constrained  approximation  to  g.  If  such  a  solution  is 
feasible  then  f  can  be  found  using  the  following  iterative  adjustment 
procedure. 

For  the  k-th  iteration  eqn.  15  is  written  as 

H  fk  *  V  <2C 


where  g  is  the  k-th  constrained  response  vector  corresponding  to 
k 

the  k-th  iteration  of  the  tap  weight  vector  f^.  For  k  *  0,  choose 
the  least-square-error  solution  f  *  f  =  S g  then  compute  gQ  from 
eqn.  20  as  gQ  *  HfQ.  The  tap  weight  vector  is  adjusted  or  "stepped" 
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using  the  relation 

fk  +  1  =  fk  +  Afk’  (21) 

k  =  0,  1,  2,  .  .  . 

where  &f^  is  given  by 

Afk  -  S  Agk  (22) 


and  where  Ag^  itself  depends  on  the  following  combination  of  the 
errors 

0  if  b(m)  £  gk(m)  —  a(™) 

-  a[g  (m)  -  a(m) ]  if  g,  (m)  >  a(m)  (23) 

k  K 

-  a[gk(m)  -  b(m) ]  if  gk(m)  <  b(m) 

The  unspecified  positive  constant  a  in  eqn.  23  must  be  sufficiently 
small  to  ensure  a  stable  iterative  process.  If  a  solution  is  feasible 
then  to  guarantee  convergence  it  has  been  proved, A  for  a  similar 
iterative  scheme,  that  Ct  must  be  equal  to  or  less  than  half  the 
reciprocal  of  the  norm  of  H . 

In  the  special  case  where  the  envelope  is  collapsed,  that  is 
a(m)  =  g(m)  =  b(m)  for  m  =  1,  2,  3,  .  .  .  M,  iteration  is  not  possible 
and  the  appropriate  solution  is  the  least-square-error  solution  f. 

If,  on  the  other  hand,  [a(m)  +  b(m)]/2  =  g(m)  and  a(m)  -  b(m)  =  2e 
for  m  =  1,  2,  3,  .  .  .M  where  E  is  the  minimax  error,  then  the 
envelope-constrained  design  should  iterate  the  tap  weights  from  their 
least-square-error  values  to  the  minimax-error  values.  In  this  sense, 
the  envelope-constrained  algorithm  relates  the  least-square  design 
to  the  minimax  design. 


V-ll 


-12- 
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5.  Summary  and  Conclusions.  When  a  finite-length  transversal  filter 

is  used  as  an  equalizer  for  a  discrete-time  system  generally  it  is  only 
possible  to  equalize  the  system  approximately.^  Three  different  4 

algorithms  have  been  presented  to  compute  the  tap  weights  of  the 

filter  subject  to  different  error  criteria  and  constraints.  Each  of 
these  design  algorithms  was  presented  in  matrix  form  suitable  for 
digital  computation.  While  the  least-square-error  design  is  fast,  it 
may  result  in  large,  unconstrained  equalization  (deconvolution)  errors. 

The  minimax  design  provides  optimum  control  of  individual  deconvolution 
errors.  The  envelope-constrained  algorithm  is  an  intermediate  design 
which  offers  greater  individual  error  control  at  the  expense  of  larger 
mean-square  error.  For  practical  equalizer  design,  it  is  probably 
worthwhile  to  evaluate  all  three  designs  and  then  select  the  one  that 
provides  the  most  suitable  approximate  equalization. 
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abstract 

The  time-domain  response  of  a  linear  system  in  cascade  with  a 
filter  is  determined  by  the  convolution  of  their  respective  impulse 
responses.  Thus,  if  the  system  to  be  equalized  is  causal  with  impulse 
response  h(t),  and  there  is  a  particular  desired  response  g(t) ,  then 
the  filter's  impulse  response  f(t)  must  satisfy, 

ft 

f (x)h(C-x)dX  =  g(t).  (1) 

o 

Since  the  unknown  filter  response  appears  under  the  integral  sign, 

(1)  is  an  integral  equation  for  f(t).  The  solution  of  (1)  requires 
that  an  inverse  convolution  or  deconvolution  be  performed  to  find 
f(t)  when  h(t)  and  g(t)  are  specified. 

Deconvolution  is  accomplished  by  direct  solution  of  the 
time-domain  convolutional  integral  equation  (1).  Previous  solutions 
to  this  integral  equation  have  been  found  by  the  method  of  colocation, 
or  by  minimizing  the  square  of  the  solution  error,  or  by  using 
envelope-constrained  procedures.  The  proposed  minirr.ax  or  Chebvshev 
solution  minimizes  the  maximum  solution  error  and,  thereby,  provides 
a  means  to  uniformly  control  deconvolution  errors. 

The  mini  max  solution  is  found  using  an  iterative  procedure 
known  as  the  second  algorithm  of  Reir.ez. 
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1.  INTRODUCTION 

1.1  Equalization 

Equalization  is  a  practical  and  important  problem  in  communica¬ 
tion  systems  engineering.  The  purpose  of  equalization  is  to  compensate 
for  the  linear  distortions  an  electrical  signal  undergoes  when  passed 
through  a  linear  system  or  communication  channel.  This  report  will 
briefly  describe  current  time-domain  equalization  techniques,  and  then 
present  a  new  approach  to  the  problem  called  minimax  time-domain 
deconvolution. 

Assume  that  the  system  or  channel  to  be  equalized  has  impulse 
response  h(t)  and  that  after  equalization  the  desired  impulse  response 
is  g(t).  Equalization  can  be  accomplished  using  a  filter  with  impulse 
response  f(t)  in  cascade  with  h(t)  as  shown  in  Fig.  1.1. 


SYSTEM  FILTER 


Figure  l.l  Equalization  Scheme 
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From  linear  system  theory  the  combined  response  g(t)  is  given 
by  the  following  convolutional  integral. 


g(t)  = 


(  h(x)f(t-x)dx  =  h(t)*f(t) 


o 


(1.1) 


The  integral  in  (1.1)  is  evaluated  between  zero  and  t  because  h(t)  and 
f(t)  are  assumed  to  be  causal. 

Mathematically,  the  continuous  functions  on  an  interval  [a,b] 
represent  an  abstract  vector  space  and  the  operator  h*(*),  representing 
convolution,  is  a  linear  transformation  on  the  space  of  continuous 
functions.  Equation  (1.1)  can  be  viewed  as  a  linear  transformation 

from  the  space  of  filter  responses  f(t)  to  the  space  of  desired 
responses  g(t)  as  shown  in  Fig.  1.2. 


Figure  1.2  Convolutional  Transformation 

Although  [a,bj  could  be  (-  °°,®),  in  practice  this  interval  is  finite. 

In  Fig.  1.2  there  is  a  particular  vector  g(t)  in  the  range 
which  is  the  most  desirable  response.  The  process  of  searching  for  the 
appropriate  vector  f(t)  in  the  domain  to  satisfy  (1.1)  is  known  as 


ppwwwp 
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One  mathematical  representation  of  the  impulse  response  for  a 
transversal  filter  is  a  series  of  N  contiguous  rectangular  pulses  of 
heights  f  and  of  width  T.  However,  an  alternate  form  which  greatly 
simplifies  the  following  discussion  is  the  Dirac  Comb, 

N 

f(t)  =  l  f  6(t-nT)  (1.2) 

n=l  n 

where  the  f  are  the  filter  tap  weights  and  6(t)  is  the  unit  impulse 
n 

function. 

Substituting  the  Dirac  Comb  representation  (1.2)  into  (1.1) 
and  integrating  yields  the  summation, 

N 

g(t)  =  l  f  h(t-nT)  (1.3) 

n=l 

for  which  h(t)  is  a  given  system  or  channel  impulse  response  and  the 

f  are  variable  scalar  coefficients.  The  deconvolution  problem  is 
n 

concerned  with  adjusting  the  N  tap  weights  f^  to  force  equality  in 

(1.3). 


1 . 3  Approximate  Solutions  to  Deconvolution  Problems 

The  linear  transformation  described  in  (1.1)  is  well  known 
in  mathematical  and  physical  sciences  as  an  integral  transformation. 
Finding  the  vector  f(t)  which  satisfies  the  transformation  for  a 
particular  g(t)  is  often  difficult  or  impossible.  The  difficulties 
are  present  because  often  the  image  space  of  the  transformation  is 
smaller  than  the  range  space.  This  means  that  as  the  transformation 
operates  on  each  vector  f(c)  in  the  domain  space,  the  set  of  vectors 
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g(t)  produced  by  the  transformation  is  only  a  subset  of  the  set  of  all 
possible  vectors  (see  Fig.  1.4). 


SPACE  OF  ALL  POSSIBLE 
DESIRED  FUNCTIONS 


Figure  1.4  Approximation  Problem  in  Deconvolution 

If  a  particular  desired  response  gQ(t)  is  not  mapped  by  the 
transformation  as  it  ranges  over  all  possible  filter  functions  (see 
Fig.  1.4)  then  it  is  desirable  to  find  a  vector  g^(t)  in  the  image 
space  which  is  in  some  way  close  to  gp(t).  Thus,  g^(t)  Is  an 
approximation  to  gQ(t)  and  there  exists  an  f(t)  in  the  domain  which 
is  mapped  to  g  (t)  by  the  transformation  h*(-)- 

A 

The  equality  stated  by  (1.3)  is  satisfied  only  if  the  choice 
of  g(t)  can  be  written  as  a  linear  combination  of  the  h(t-nT). 
Deconvolving  (1.3)  to  determine  the  values  of  the  N  tap  weights  f^ 
is  straightforward  only  when  g(t)  lies  in  the  span  of  the  N  time 
shifted  functions  h(t-nT).  This  report  will  investigate  possible 
approximate  solutions  for  the  more  general  case  when  g(t)  is  not  . 
spanned  by  the  right-hand  side  of  (1.3). 
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1.4  Approximation  Theory 

Referring  to  Fig.  1.4  it  is  desired  to  pick  an  approximate  vector 
g^(t)  so  that  the  distance  between  it  and  the  desired  vector  gQ(t)  is 
small.  To  determine  a  best  approximation  it  is  first  necessary  to  define 
a  measure  of  closeness  between  vectors  of  an  abstract  vector  space.  For 
example,  it  is  not  clear  what  the  distance  between  the  continuous 
functions  cos(t)  and  sin(t)  is  on  the  interval  [0,  2tt ] . 

Intuitively,  a  distance  operator  d(.,0  on  pairs  of  vectors 
must  have  certain  properties.  The  operator  assigns  a  real  number  to  the 
vector  pair  such  that  the  distance  between  the  two  vectors  is  greater 
than  zero  unless  the  two  vectors  are  equal.  If  the  vectors  are  equal, 
then  the  operator  assigns  zero  to  be  the  distance.  Furthermore,  the 
distance  operator  must  be  independent  of  direction  and  obey  the  triangle 
inequality. 

Summarizing  for  the  vectors  x,  y,  z 

d(x,y)  _>  0  ,  x  /  y 

d(x,x)  =  0 

(1.4) 

d(x,y)  =  d (y , x) 
d(x,z)  <_  d(x,y)  +  d(y,z) 

A  well-known  distance  operator  on  the  real  numbers  is  the 
absolute  value  function  which  obeys  the  properties  of  (1.4)  for  any  real 
numbers  x,y,z.  For  abstract  vector  spaces  a  general  class  of  operators 
called  norms  have  the  properties  of  a  distance  operator.  A  norm  assigns 
a  real  number  to  a  vector  and  is  denoted  by  ]|*i|. 


The  norm  has  the 


ft 
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following  properties  for  any  vectors  x,  y 

I  M  I  =  0  ,  x  =  0 

I |x||  >  0  ,x^0 

I laxl I  =  I a I  *  1 1 x! I  .  a 
I |x+y| !  1  I lx! I  +  I |y| I 


It  can  easily  be  shown  that  | |x-y| |  is  a  distance  operator  on  the  pair 
of  vectors  x,  y. 

Referring  to  Fig.  1.4  and  using  the  norm  as  the  measure  of 
distance,  it  is  desired  to  minimize  the  distance  between  the  desired 
response  g  (t)  and  the  approximate  response  g  (t).  For  an  interval 

A 

[a,b]  this  procedure  is  represented 

min  | |gD(t)  -  gA(t) | ]  (1.6) 

and  is  determined  by  examining  (1.6)  for  all  possible  g^(t). 

Since  gA(t)  can  be  represented  as  a  linear  combination  of  the  N 
functions  h(t-nT)  (1.6)  becomes 

N 

min  |]gD(t)  -  l  f ^h ( t-nT) | | .  (1.7) 

n=l 


It  is  convenient  to  represent  the  arguments  of  (1.6)  and  (1.7)  as  an 
error  function 
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N 

r (t)  -  g(t)  -  l  f  h(t-nT)  (1.8) 

n=l 

where  a  best  approximation  is  recognized  when  the  norm  of  r(t)  is 
minimized.  An  important  point  is  that  r(t)  is  a  function  of  the  N  tap 
weights  f^  and  that  |  |r(t) j  |  will  vary  as  f  vary. 

A  norm  that  is  widely  used  as  a  distance  measure  in  the 

engineering  sciences  is  the  Euclidean  norm.  The  Euclidean  norm  of  a 

2  2 

x  +  y  and  for  two  vectors, 

v^  =  (x2»y-]^  and  v2  =  (x2’y2^  the  norm  of  the  difference  is 

I!vi-v2II  =/(xi~x2)2  +  {yi~y2)2  (1‘9) 


Equation  (1.9)  can  be  generalized  to  an  n-dimensional  space  Rn  . 

For  a  vector  x  =  (x, , . . . ,x  ) ,  x  ERn 
i  n 


(1.10) 


There  is  an  analogous  representation  for  the  Euclidean  norm  operating 
on  the  abstract  vector  space  of  continuous  functions  on  an  interval 
[a,b].  For  ^(t)  continuous  on  [a,b] 


ll*(t)|| 


and  for  <t>^(t)  and  d>2 ( c )  continuous  on  [a,b] 


(1.11) 
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11^(0  -  4>2 ( t)  1 1 


(1.12) 


In  electrical  engineering  (1.11)  is  the  RMS  (root  mean  square)  value  of 
<J>(t). 

Using  the  Euclidean  norm  as  defined  in  (1.12)  it  is  now 
possible  to  calculate  a  distance  between  sin(t)  and  cos(t),  the  functions 
of  the  previous  example, 


| | sin  t  -  cos  t| |  = 


cost)  dt  =  1. 


(1.13) 


Note  that  the  value  of  the  norm  in  (1.13)  is  related  to  the  area  between 
the  curves  sin(t)  and  cos(t)  in  the  interval  [0,2^].  This  is  an 
important  characteristic  of  the  Euclidean  norm.  Minimizing  the 
Euclidean  norm  is  in  essence  minimizing  the  area  between  an  approxima¬ 
tion  curve  and  the  desired  curve. 

Minimization  of  the  Euclidean  norm  is  the  approximation  tech¬ 
nique  known  as  least  squares.  The  characteristic  feature  of  least 
squares  solutions  is  that  the  area  between  the  desired  curve  and  the 
approximate  curve  is  minimized.  A  significant  drawback  of  the  method 
is  that  individual  errors  are  not  bounded.  Often  a  least  squares 
solution  has  large  deviations  from  the  desired  curve  at  certain  points. 
For  example,  the  Gibb's  phenomenon  observed  with  finite  Fourier  series 
solutions  produces  sharp  spikes  in  the  error  curve  at  discontinuities 


of  the  desired  curve. 


One  case  where  large  deviations  in  the  error  curve  are  particu¬ 
larly  undesirable  is  in  digital  signal  detection.  Since  the  detection 
scheme  is  often  a  simple  level  detector  the  presence  of  large  deviations 
in  the  approximation  curve  could  result  in  erroneous  detection  of 
information  during  transmission. 

An  approach  to  limiting  unbounded  errors  in  the  least  squares 
approximation  scheme  has  been  explored  in  the  literature.  The  method 
is  called  envelope-constrained  deconvolution  and  it  is  an  iterative 
technique  which  is  initiated  by  computing  the  least  squares  approximation. 
In  successive  iterations,  the  largest  errors  are  reduced,  if  possible, 
until  they  are  bounded  by  an  envelope  defined  around  the  desired  curve 
as  in  Fig.  1.5. 


g(t) 


Figure  1.5  Envelope  Constrained  Deconvolution 

Pushing  errors  down  in  one  place  may  push  errors  up  in  other  places. 
Generally,  mean  square  error  is  greater  when  envelope  constraints 
are  imposed  to  suppress  large,  individual  errors.  The  envelope-, 
constrained  approximation  is  a  practical  approach  to  restricting  the 


magnitude  of  the  error.  The  difficulty  is  that  it  is  not  generally  known 
how  narrow  the  envelope  can  be  made,  and  there  is  an  important  question 
of  solution  feasibility. 

1.5  Minimax  Approximation 

Here,  a  different  approach  to  the  deconvolution  problem  is 
proposed,  namely,  the  minimax  solution.  The  minimax  solution  is  a  best 
bounded  approximation  and  is  derived  as  a  characteristic  of  using  the 
Chebyshev  norm  as  a  distance  measure  instead  of  the  Euclidean  norm. 

The  Chebyshev  norm  of  a  vector  x  =  (x, ,...,x  )  which  is  a 

1  n 

member  of  Rn  is 


max 

ljfUn 


(1.14) 


The  analogous  representation  of  (1.14)  for  continuous  functions  on  an 
interval  [a,b]  is 


I  !<Kt)  |  |  =  a™*b  |<KO|  (1.15) 

where  <p(t)  is  a  continuous  function  of  [a,b].  The  minimax  approximation 
comes  from  minimizing  (1.14)  or  (1.15). 

The  remainder  of  this  report  deals  with  further  discussion  of 


the  minirax  approximation  and  its  applications  to  the  design  of  trans¬ 
versal  filter  equalizers.  Future  reference  to  the  norm  or  the  use  of 
the  norm  will  imply  the  Chebyshev  norm  unless  otherwise  noted. 


24. 


(2.2) 


The  matrix  of  coefficients  A  =  [a  ]  on  the  left-hand  side  of 

mn 

(2.2)  is  a  linear  transformation  on  the  N-ditnensional  vectors 
x  =  (x^,  ...,  x^)  ,  xCRn  to  the  M-dimensional  vectors  y  =  (y^,  ...»  v^) , 
ye  Rm.  The  vector  b  =  (b^,  •••t  b^)  on  the  right-hand  side  of  (2.2)  is 
a  constant  vector  for  which  it  is  desired  to  find  the  vector  x  which 
under  the  transformation  A  maps  to  b: 

Ax  =  b.  (2.3) 

Equation  (2.3)  is  equivalent  to  (2.2). 
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The  existence  and  the  number  of  solutions  to  (2.3)  depends  upon 
the  relationship  between  the  number  of  equations  and  the  number  of 
unknowns.  Three  possibilities  exist,  they  are  M<N,  M=N,  and  M>N  where  M 
is  the  number  of  equations  and  N  is  the  number  of  unknowns. 

For  each  case,  the  transformation  A  is  mapping  vectors  from 
N-dimensional  space  Rn  to  M-dimensional  space  R31.  In  the  case  M<N  the 
transformation  maps  from  a  larger  space  to  a  smaller  space  and  if  the  rows 
of  the  matrix  are  linearly  independent,  then  for  any  particular  vector  b 
which  belongs  to  Rm  there  are  always  infinitely  many  vectors  xeRn  such 
that  b  is  mapped  by  these  vectors  under  the  transformation.  Referring  to 
Fig.  1.2  in  Chapter  1,  it  is  intuitively  easy  to  see  how  a  domain  space 
of  larger  dimension  would  fill  up  many  times  over  a  range  space  of 
smaller  dimension. 

In  the  case  M=N  the  transformation  maps  from  a  domain  space  of 
equal  dimension  to  the  range  space.  From  linear  algebra  it  is  known 
that  a  system  of  linear  equations  with  the  same  number  of  equations  as 

k 

unknowns  has  a  unique  solution  if  the  equations  are  linearly  independent. 

Finally,  for  the  case  M>N  the  transformation  operates  on  a 
domain  space  which  is  smaller  than  the  range  space.  As  in  Fig.  1.4  of 
Chapter  1  it  can  be  seen  that  the  transformation  of  the  domain  space  does 
not  fill  up  the  range  space.  That  is,  the  transformation  of  the  domain 
is  a  subspace  of  the  range.  If  the  vector  b  in  (2.1)  -  (2.3)  lies  in 
the  image  space  of  the  transformation  A,  then  a  solution  exists. 
Alternatively,  if  b  docs  not  lie  in  the  image  space  of  A  then  no  exact 
solution  exists.  When  an  exact  solution  does  not  exist  the  best' 

* 

0  nan,  H. ,  Linear  Algebra,  Harcourt  Brace  Jovanovich,  1976 
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solution  to  (2.1)  -  (2.3)  is  found  by  selecting  a  vector  which  lies  in 
the  image  space  of  A  and  is  a  closest  approximation  to  b. 

When  the  number  of  equations  in  a  linear  system  exceeds  the  number 
of  unknowns  the  system  of  equations  is  called  overspecified.  To  find  the 
best  approximation  b  to  the  vector  b  the  Chebyshev  norm  will  be  used  as  a 
distance  measure.  If  £  is  the  smallest  distance  that  can  be  achieved 
between  the  two  vectors,  then 

e  -  | |b  -  b |  |  (2.4) 

Since  b  is  mapped  by  the  transformation  A  for  some  choice  of  x  (2.4) 
becomes 

£  =  min  | J  b  -  Ax  | |  (2.5) 


2.2  Solving  Two  Equations  in  One  Unknown 

Often  only  approximate  solutions  exist  for  overspecified 
systems  of  linear  equations.  This  section  and  the  rest  of  this  chapter 
will  be  concerned  with  describing  the  characteristics  of  properties  and 
the  minimax  approximation. 

The  linear  equation  2x  =  1  has  an  exact  algebraic  solution, 
namely,  x  =  1/2.  The  equation  2x  =  1  is,  of  course,  the  most 
elementary  case  of  a  linear  system  of  equations.  The  solution  to  the 
equation  2x  =  1  can  be  found  graphically  by  defining  an  error  equation 
r1(x),  such  that  r^(x)  =  |2x  -  l|.  Plotting  r^(x)  in  Fig.  2.1  shows 
that  the  solution  to  the  equation  2x  =  1  is  found  where  r^(x)  intersects 
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the  x-axis.  That  is,  when  the  solution  has  been  found  the  error  is  zero. 
In  the  definition  of  r^(x)  magnitude  brackets  are  used  since  it  is  only 
the  magnitude  of  the  error  which  is  of  concern. 


r  (x) 


Figure  2 . 1  Solution  of  a  Linear  Equation  in  One  Unknown 

If  a  system  of  linear  equations  has  more  linearly  independent 
equations  than  unknowns,  then  an  exact  solution  does  not  exist.  This 
is  so,  for  example,  with  the  system: 

2x  =  1 

(l/3)x  =  1  (2.6) 

Each  separate  equation  has  a  unique,  exact  solution  (1/2  and  3 
respectively)  which  can  be  found  either  algebraically  or  by  the 
graphical  method  of  Fig.  2.1.  However,  there  is  no  solution  which 
simultaneously  satisfies  both  equations  in  (2.6)  exactly. 


Figure  2.2  which  is  a  plot  of  the  error  equations  r^(x)  and 
^(x)  associated  with  (2.6).  An  intuitive  approximate  solution  to  (2.6) 
would  be  x  =  6/7. 


!  r^(x)  =  |2x-l( 


r„(x)  =  | (1/3) x— 1 [ 


6/7  1 


Figure  2.2  Approximation  to  Two  Linear  Equations  in  One  Unknown 


Plugging  the  value  x  =  6/7  into  the  left-hand  side  of  (2.6)  yields 


2x  =  2(6/7)  =  ly 


=  j  (6/7)  =  2/7 


(2.7) 


Note  that  the  proposed  solution  x  =  6/7  produces  values  in  (2.7) 
which  are  in  error  with  respect  to  the  right-hand  side  of  (2.6)  by 
5/7  for  each  equation.  Furthermore,  the  value  x  =  6/7  is  where  the 
functions  r^(x)  and  t^Cx)  intersect  in  Fig.  2.2  r.d  r^(x) 

also  intersect  at  x=0,  however,  the  error  produced  for  each  equation 
when  zero  is  used  in  the  left-hand  side  of  (2.7)  is  r^  =  r^  =  1  which 
is  greater  than  the  error  produced  by  x  =  6/7). 
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Now  the  minimax  approximation  to  (2.6)  will  be  computed.  The 
system  of  equations  in  (2.6)  are  of  the  form  of  (2.3)  where  the  matrix 


*  2  ' 

~1~ 

of  coefficients  is  A  = 

.1/3. 

and  the  vector  b  = 

1 

Using  the  norm 


as  the  distance  measure  between  the  approximation  Ax  and  the  exact 
solution  b  yields 


~l' 

'  2x 

"l  -  2x~ 

1 

II  =  II 

,  1 

1 

L_  — 

L  3XJ 

_L  '  3XJ 

(2.8) 


The  minimax  solution  is  found  by  minimizing  the  norm  of  (2.8).  Using 
the  Chebyshev  norm  and  denoting  C  as  the  minimax  error, 


~1 

-  2x 

“  ,  “I 

£  =  min  |  | 

1 

1 

-  :rx 

|  |  =  min 

max  {  |  l-2x  |  ,  |l-^x|  } 

3  J 

r 


=  min  jjaax  (r^(x),  r^x)} 


(2.9) 


The  function  max  (r^(x),  ^(x)}  in  (2.9)  is  the  crosshatched  portion 
of  the  graph  in  Fig.  2.2,  and  the  minimum  value  of  this  function  is 
found  when  x  =  6/7  as  shown.  The  value  of  the  error  when  x  =  6/7  in 
(2.9)  is  e  =  5/7. 

The  intuitive  approximation  x  =  6/7  was  chosen  because  it  was 
the  value  which  minimized  the  error  for  hoth  equations  in  (2.6).  Note 
that  moving  either  to  the  left  or  right  of  x  =  6/7  in  Fig.  2.2  increases 
the  error  for  one  of  the  two  error  equations.  The  minimax  approximation 
is  equivalent  to  the  intuitive  approximation,  however,  the  procedure  for 
finding  the  minimax  approximation  can  be  formulated  graphically. 
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2.3  Solving  M  Equations  in  One  Unknown 

The  graphical  technique  for  finding  the  minimax  approximation 
to  a  set  of  two  equations  in  one  unknown  can  be  expanded  to  a  larger  set 
of  equations.  Consider,  for  example,  the  set  of  equations 

2x  =  1 
l/3x  =  1 
5/2x  =  5 

5/8x  =  5/2  (2.10) 

The  associated  error  equations  are 

r1(x)  =  |2x  -  l| 

r2(x)  *  1-5*  "  x! 

r3(x)  =  |-|x  -  5 | 

r4(x)  =  ||x  -  ||  (2.11) 

Defining  an  error  vector  r(x)  which  has  the  four  components 
r^(x),  r2(x),  r3(x),  r^(x)  it  is  possible  to  find  the  norm  of  r(x)  from 
the  graph  of  the  error  functions  r^,  r^,  r^.  The  norm  of  r(x)  is 

l  i  r  ^  x )  |  j  max 

ind  it  is  the  crosshatched  region  of  Fig.  2.3.  From  Fig.  2.3  it  is  easy 
to  pick  out  the  minimum  value  of  (2.12)  and  the  value  of  x  associated 
to  it.  These  values  are  min  ||r(x)||  =  5/3  and  x  =  4/3.  Thus  x  3  4/3 
:s  t.he  minimax  approximation  to  (2.10). 
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Several  characteristics  of  the  graphical  solution  to  (2.10)  are 
important  and  deserve  further  discussion.  First,  the  minimax  solution 
x  =  4/3  is  found  at  the  intersection  of  two  of  the  error  equations  in 
(2.11).  Referring  to  Fig.  2.3  the  particular  equations  are  r^(x)  = 

1 2x-l  |  and  r^(x)  =  |-|  x  -  -|j  . 


Figure  2.3  Approximation  to  Four  Linear  Equations  in  One  Unknown 

Note  that  all  other  error  equations  lie  below  the  intersection  of  r^(x) 
and  r^(x)  at  x  =  4/3.  Thus,  the  maximum  error  has  been  minimized  for 
all  four  equations.  If  an  approximation  is  chosen  slightly  to  the  left 
or  right  of  x  =  4/3  then  either  the  value  of  r^  increases  or  the  value 
of  r^  increases.  For  x  =  4/3  the  error  for  each  equation  has  been 
bounded  and  minimized. 


32. 


Second,  the  solution  occurs  at  the  intersection  of  two  lines  which 
are  of  opposite  slope  (see  Fig.  2.3).  This  is  a  very  important  property 
and  is  related  to  the  fact  that  the  Chebyshev  norm  of  r(x)  is  a  convex 
function. 

A  convex  function  is  defined  as  being  a  function  cj) (x)  such  that 
for  any  three  values  of  x,  a<b<c,  <{>(b)  04> (a)  +  (1-0)  (c)  where  0  varies 

from  zero  to  one.  Figure  2.4  shows  a  convex  function. 


A  convex  function  always  has  a  global  minimum.  Since  the  norm  of  r(x) 
is  a  convex  function  (see  Fig.  2.3)  a  global  minimum  can  be  found  for 
it,  and  the  minimax  approximation  is  guaranteed  to  exist. 

A  change  in  the  sign  of  the  slope  of  ||r(x)||  at  the  minimax 
solution  is  an  indication  that  a  global  minimum  of  the  convex  function 
J | r (x) | |  has  been  found.  This  can  be  seen  by  moving  along  ||r(x)||  in 
Fig.  2.3.  If  a  point  of  intersection  is  made^of  two  lines  with  the  same 
sign  in  slope,  it  indicates  that  ||r(x)||  can  be  further  decreased. 
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2.4  Characterization  of  the  Minimax  Approximation 

The  two  properties  described  at  the  end  of  section  2.3  apply  to 

overspecified  systems  of  linear  equations  with  more  than  one  unknown. 

The  generalization  of  these  two  properties  specifies  the  minimax  solution 

* 

and  is  known  as  the  Characterization  Theorem.  In  section  2.3  which 
considered  systems  of  linear  equations  in  one  unknown  the  minimax  approx¬ 
imations  were  found  at  the  intersection  of  two  of  the  equations. 
Furthermore,  the  slopes  of  the  two  equations  were  of  opposite  sign. 

In  general,  a  system  of  M  linear  equations  in  N  unknowns  with 
M>N,  the  minimax  approximation  is  found  at  the  intersection  of  some  subset 
of  N-f-1  of  the  M  equations.  Generalizing  to  N  unknowns  the  property  that 
the  minimax  approximation  is  found  at  the  intersection  of  two  equations 
which  have  slopes  of  opposite  sign  requires  the  introduction  of  the 
concepts  of  a  convex  set  and  a  convex  hull 

A  convex  set  is  a  set  of  vectors  such  that  for  any  two  vectors 
v  and  w  contained  in  the  set,  the  set  also  co  .tains  all  vectors  of  the 
form  9v  +  (l-G)w  where  0  is  a  scalar  which  varies  from  zero  to  one. 

2 

Figure  2.5  shows  some  examples  of  convex  and  non-convex  sets  in  R  . 


CON\' EX  CONVEX  NOT  CONVEX 


Figure  2.5  Convex  Sets 


* 

Cheney,  E.W.,  Introduction  to  Approximation  Theory,  McGraw  Hill, 
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The  convex  hull  of  a  set  S  which  itself  is  not  convex  is  the  smallest 
set  which  contains  S  and  is  convex. 


Now  it  is  possible  to  generalize  the  fact  that  for  systems  with 

one  unknown  the  slopes  are  of  opposite  sign.  If  the  row  of  the  matrix 

on  the  left-hand  side  of  (2.2)  is  denoted  A  then  to  determine  whether  a 

m 

particular  intersection  of  N+l  error  equations  determines  the  minimax 
approximation  depends  upon  a  property  of  the  N+l  associated  rows  of  the 


matrix  A  (i.e.,  associated  by  the  expression  r  (x)  =  lb  -I  a  x  I 

m  '  m  ,  mn  n' 

n=l 

where  A  =  (a  . , . . . , a  ..) ) .  Denoting  these  N+l  rows  as  A  A 

m  ml  mN  m(l)  m(N+l) 

the  condition  states  that  the  minimax  approximation  exists  when  zero 
belongs  to  the  convex  hull  of  the  set 


{°m(l)Am(l) ,am(N+l)Am(N+l) }  * 


(2.12) 


This  is  where  a  .a  are  the  signs  +  1  determined  by  the 

m(l)  m(N+l)  &  * 

equation 


N 

a  =  sgn  (b  -  7  a  x  } 
m  °  m  ^  mn  n 
n=l 


(2.13) 


f0r  Am  =  ^ml---**^- 

An  equivalent  statement  of  the  condition  that  zero  life  in  the 
convex  hull  of  (2.12)  is  that 


N+l 


0  1 


n=l 


m(n)  m(n)  in(n) 


See  section  A. 2. 


(2.14) 
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where  the  0  .  .  are  scalar  constants  and 
m(n) 


N+l 

0  .  .  >  0  and  J  0  .  .  =  1 
m(n)  —  m(n) 

n=i 


(2.15) 


In  the  case  of  one  unknown  the  vectors  A  are  one-dimensional 

m 

and  thus  belong  to  the  real  line  R .  Two  of  the  vectors  determine  a  mini¬ 
max  approximation  only  if  their  slopes  are  of  opposite  sign.  For  example, 
the  two  vectors  could  be  -5  and  3  which  belong  to  R .  From  the  definition 
of  the  convex  hull,  it  is  true  that  zero  is  contained  in  the  convex  hull 
of  the  set  (-5,  3}  . 

A  proof  of  the  facts  that  the  minimax  approximation  is  found  in 

some  subset  of  N+l  of  the  M  error  equations  and  that  zero  is  contained 

in  the  convex  hull  of  the  vectors  A  ,  .  associated  with  the  subset  of 

m(n) 

N+l  error  equations  can  be  found  in  the  appendix  in  sections  A.l  and  A. 3. 


3.  USING  A  TRANSVERSAL  FILTER  AS  A  DISCRETE  TIME-DOMAIN  EQUALIZER 


3.1  Convolution  and  the  Transversal  Filter 


As  discussed  in  section  1.2  a  transversal  filter  has  an  impulse 
response  which  can  be  represented  as  in  (1.2) 


N 

f(t)  =  J  f  S(t-nT)  (3.1) 

n=l  n 


Since  the  impulse  response  of  the  filter  is  finite  in  length,  it  is 
known  as  a  Finite  Duration  Impulse  Response  (FIR)  filter.  Referring 
to  Fig.  3.1  a  transversal  filter  is  composed  of  N  sections  each  of  which 
consists  of  a  delay  T  and  a  scalar  multiplier  f^.  An  incoming  signal 
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is  first  delayed  and  then  multiplied  by  the  scalar  weight  of  the  first 
section.  The  next  section  further  delays  the  unweighted  signal  and  then 
multiplies  it  by  its  own  scalar  weight.  The  signal  traverses  all  N  sections 
of  the  filter  and  the  output  of  the  filter  is  a  summation  of  the  N  delayed 
and  weighted  versions  of  the  input  signal.  The  scalar  weights  are  variable 
and  can  be  manipulated  either  electrically  or  mechanically. 


Figure  3.1  Transversal  Filter 

If  the  Input  to  a  transversal  filter  is  a  discrete-time  signal 
x(t)  vftjach  only  takes  non-zero  values  at  equally  spaced  discrete  points 
of  time,  then  x(t)  can  be  represented  by  x(nT)  where  T  is  the  time 
between  non-zero  values  of  x(t)  and  n  is  an  integer.  The  transversal 

. . .  .  ..  ua  <x  aiscre . 

signal  with  non-zero  values  f(nT). 

The  output  of  a  transversal  filter  with  input  x(nT)  is  the 
convolution  sum  of  x(nT)  and  f(nT).  Calling  y(mT)  the  output  of  the 
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filter  for  m  an  integer  and  assuming  that  the  value  of  T  is  one  yields 

N 

y(m)  =  l  x(m-n)f (n) .  (3.2) 

n=l 

The  convolution  sum  of  (3.2)  is  from  n=l  to  n=N  since  f(n)  is  of  finite 
length  N. 

For  a  discrete-time  system  or  communication  channel  which  has 
impulse  response  h(nT)  =  h(n)  for  T  =  1  the  combined  impulse  response  of 
h(n)  in  cascade  with  f(n)  is  the  convolution  sum  of  the  two.  If  h(n) 
is  causal  starting  at  n=l  and  of  finite  length  L,  then  the  combined 
response  g(m)  is 


N 

g(m)  =  l  h(m-n+l)f (n)  m=l,...,M  (3.3) 

n=l 


where  the  length  of  g(m)  is  M  =  N+L+l.  The  operation  of  convolution  sum 
as  in  (3.2)  and  (3.3)  is  often  represented  by  an  asterisk  so  that  (3.3) 
becomes 


g(n)  =  h(n)*f(n).  (3.4) 

3.2  Solving  the  Equalization  Problem 

The  purpose  of  the  discussion  to  this  point  has  been  to  explain 
how  the  transversal  filter  can  act  as  an  equalizer  for  a  system  or  a 
communication  channel.  A  system  or  channel  may  have  an  undesirable 
impulse  response  because  in  some  way  it  introduces  linear  distortion 
(i.e.,  alters  the  waveshape)  to  signals  which  traverse  it.  Equalization 


is  the  process  of  changing  the  impulse  response  of  the  system  or  channel 
through  the  addition  of  a  transversal  filter  as  an  equalizer.  Thus,  the 
combination  of  the  system  or  channel  with  the  filter  has  a  more  desirable 
reponse. 

If  the  system  or  channel  response  is  measured  and  called  h(n) 
and  the  desired  response  is  specified  as  g(m)  then  the  problem  is  to 
determine  the  f(n)‘s  which  satisfy  (3.3).  Equation  (3.3)  is  of  the  form 
(2.1)  in  section  2.1  and  can  be  written  in  the  form  of  (2.2) 


Ml) 

h(2) 


h(N) 


0 

h(l) 


0 

0 


h(N-l)  h(N-2) 


h(M-N)  h(M-N-l)  h(M-N-2) 
h (M-N+l)  h(M-N)  h(M-N-l) 
0  h (M-N+l)  h(M-N) 

0  0  h (M-N+l) 


0 

0 


h(l) 


h(M-2N) 
h(M-2N+l) 
h(M-2N+2) 
h  (M-2N+3), 


h (M-N+l) 


—  — 

-  — 

f(l) 

8(1) ' 

f(2) 

8(2) 

f(N) 

g(M) 

_  — 

L-  _ 

(3.5) 


If  M  is  greater  than  N,  which  must  be  the  case  for  L>1  ,  then  (3.5)  is 
an  overspecified  system  of  linear  equations  where  the  unknowns  are  the 


From  the  discussion  in  section  2.1  it  is  known  that  an  over- 
specified  system  of  linear  equations  generally  has  no  exact  solution. 
Thus,  to  solve  the  equalization  problem  it  is  desired  to  find  a  best 
approximate  solution  to  (3.5)  which  may  be  accomplished  as  described  in 
Chapter  2. 
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4.  DISCRETE  TIME  ALGORITHM 

4.1  Ascent  Algorithm  for  One  Unknown 

In  the  equalization  problem  a  system  impulse  response  h(n)  is 
given  and  a  desired  response  g(n)  is  specified.  The  tap  weights  f(n)  are 
unknown.  As  pointed  out  in  Chapter  3  the  equalization  problem  generates 
an  overspecified  system  of  linear  equations  in  N  unknown,  the  f(n)  being 
those  unknowns.  The  algorithm  desire’  must  seek  the  minimax  approxima¬ 
tion  to  the  f(n). 

The  graphical  technique  of  Chapter  2  for  finding  the  minimax 
approximation  to  an  overspecified  system  of  linear  equations  in  one 
unknown  is  not  useful  as  an  algorithm  when  the  system  has  more  than  one 
unknown.  A  generalized  algorithm  is  needed  which  has  the  capacity  to  find 
the  minimax  approximation  to  an  overspecified  system  of  linear  equations 
in  N  unknowns. 

This  chapter  gives  an  intuitively  plausible  algorithm  for  the 

case  of  an  overspecified  system  in  one  unknown.  Then  the  algorithm  is 

generalized  to  the  case  of  N  unknowns.  The  algorithm  to  be  described 

k 

is  known  as  the  ascent  algorithm.  For  the  case  of  one  unknown  a 
graphical  description  of  the  algorithm  is  based  on  properties  of  the 
minimax  approximation  as  described  in  the  characterization  theorem  in 
section  2.4. 

Referring  to  Fig.  4.1,  a  particular  iteration  of  the  algorithm 
has  located  a  point  xq  corresponding  to  the  intersection  of  two  lines. 

Each  line  represents  an  error  equation  (see  section  2.3),  and  these  lines 


* 

Cheney,  E.W. ,  Introduction  to  Approximation  Theory, 
>v.Graw  Hill,  1966. 
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have  slopes  of  opposite  sign  satisfying  the  requirements  of  the  char¬ 
acterization  theorem.  The  characterization  theorem  states  that  zero  must 
belong  to  the  convex  hull  of  the  vectors  which  represent  the  slopes  of 
the  lines.  This  condition  is  satisfied  in  the  case  of  one  unknown  only 
when  the  slopes  of  the  two  lines  are  of  opposite  sign. 

If  xq  corresponded  to  the  global  minimax  solution  then  at  the 

point  xq  there  would  be  no  lines  above  the  two  lines  whose  intersection 

defines  x  and  the  algorithm  would  terminate.  This  is  true  because  the 
o 

minimax  approximation  is  found  at  the  intersection  of  two  lines  which 

belong  to  the  function  ||r(x)jj  =  (r^Cx)}  (see  (2.12)).  As  in 

Fig.  4.1,  if  it  is  not  the  case  that  x  is  the  global  minimax  approximation 

o 


then  a  new  line  is  chosen  such  that  the  magnitude  of  'he  chosen  line 

evaluated  at  x  is  equal  to  |]r(x  )|j  .  Now  the  aigo.it.im  .nis  cnree  iines. 
o  o 

Out  of  the  three  lines  the  algorithm  will  choose  two  and  delete  the  third. 
The  algorithm  retains  the  new  line  and  one  of  the  old  lines  which  has 
slope  of  opposite  sign  to  the  new  line. 
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Where  the  two  new  lines  intersect  is  the  new  point  of  the  algorithm. 

In  Fig.  4.1  this  point  is  x^.  The  discovery  of  defines  the  next  iteration 

of  the  algorithm  and  the  algorithm  proceeds  with  x,  as  it  did  with  x  . 

1  o 

4.2  Generalization  of  the  Ascent  Algorithm  to  N  Unknowns 


There  are  three  aspects  of  the  ascent  algorithm  which  require  gen¬ 
eralization  to  N  unknowns.  They  are  first,  how  to  ensure  that  the  convex 
hull  of  the  N+l  error  equations  contains  zero;  second,  how  to  arrive  at 
the  intersection  of  N+l  error  equations;  third,  how  to  exchange  one  of 
the  N+l  error  equations  for  the  particular  equation  that  belongs  to  the 
function  | |r (x) j | .  The  algorithm  should  solve  all  these  problems  effi¬ 
ciently  and  should  be  able  to  determine  when  the  solution  has  been  found. 


4.2.1  Forcing  zero  to  belong  to  the  convex  hull 

The  overspecified  system  of  linear  equations  representing  the 


equalization  problem  is  of  the  form  (2.2) 


311  ai2 


a21  a22  '  •  *  a2N 


(4.1) 


*  ,  T 

i’ll  MZ 


J  l_  J 


for  M>N.  For  each  row  of  the  matrix  in  (4.1)  there  is  an  associated 


error  equation 


M  -lb-  l 


m  ,  mn  n 

n=l 


. . . 
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The  algorithm  begins  by  randomly  choosing  a  subset  of  N+l  of  the  M  error 

equations.  To  denote  which  subset  of  N+l  equations  are  being  considered, 

the  error  equations  will  be  subscripted  r  /1N(x),  r  ,„.(x),  .  . 

mil;  nHz; 

rm(N+l) (x) * 


It  is  desired  that  zero  belong  to  the  convex  hull  of  the  set 

{o  ...A  ,..,a  ...A  a  A  .} where  A  is  the  mth  row  vector 

m(l)  m(l)  m(2)  m(2)  m(N+l)  m(N+l)  m 

of  the  matrix  in  (4.1)  and  a  are  the  signs  +  1  defined  by  the  equation 

m  ~ 


a 

m 


N 

sgn  {b  -Fa  x  }. 
m  ,  mn  n 
n=l 


(4.3) 


An  intuitive  understanding  of  the  purpose  of  the  can  be  obtained  by 

referring  to  the  case  of  an  overspecified  system  in  one  unknown  in 

Fig.  2.3.  Each  error  equation  in  Fig.  2.3  has  two  parts,  one  part  is  of 

negative  slope  and  the  other  part  is  of  positive  slope.  It  is  important 

to  distinguish  between  these  two  portions  of  r  (x) .  For  example,  the 
“  m 

intersection  of  the  lines  r^(x)  and  r^(x)  in  Fig.  2.3  occurs  in  two 

places  x  =  0  and  x  =  6/7.  At  x  *  0  the  two  lines  have  slopes  of  similar 

sign  and  at  x  5  6/7  the  lines  have  slopes  of  opposite  sign.  Thus,  one 

of  the  intersections  constitutes  a  miniraax  approximation  and  the  other 

does  not.  The  O  are  the  devices  which  keep  track  of  which  of  the  two 
m 

portions  of  an  error  equation  the  algorithm  is  using. 

It  is  now  possible  to  determine  what  values  of  the  O  will 

m 

make  zero  belong  to  the  convex  hull  of  the  set  (o  ,,.A  O 

mil;  mil; 

A  /,,,  i , } .  By  (2.14)  and  (2.15)  zero  will  belong  to  the  convex  hull  of 
this  set  when 

N+l 

0  -  I  0  .  .o  ,  .A  ,  . 

m(n)  m(n)  m(n) 


(4.4) 


'«  • 


0  ,  .>0  and  J  0  -  .  =  1 

min  I —  mini 


(4.5) 


By  choosing 


a  .  .  =  sgn  (0  ,  ,} 
m(n)  m(n) 


(4.6) 


and  recognizing  that  0  ,  .ct  .  .  for  0  „  s>0  is  idential  to  0  ,  .  (for 
°  m(n)  m(n)  m(n)—  m(n) 

any  positive  or  negative  value  of  0  .  .),  (4.4)  can  be  rewritten  as 

m(n) 


^  am(n)Am(n)  ^m(n)am(n)^  ^m(n) Am(n) ^ ' 


(4.7) 


Recalling  that  A  =  (a  , ,a  ~ , . . . ,a  „)  and  that  the  sum  of  the 
m  mi  mz  mN 

0 ’/  ■.  must  equal  1  allows  (4.7)  to  be  rewritten  in  the  form  of  a  matrix 
m(n) 

transformation. 


am(l)l  am(2)l  ’ 


m(N+l)l 


(4.8) 


am(l)N  am(2)N  * 


3m(N+l)N  ®m(N+l)| 


A  solution  for  the  0  .  .  exists  if  the  square  matrix  in  (4.8)  has  an 

nun) 

inverse.  The  inverse  of  the  matrix  in  (4.8)  is  defined  as  the  matrix 
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i 

♦ 

j 


such  that 


am(l)l  ata(2)l 


m(N+l)l 


am(l)N  am(2)N  .  .  .  am(N+l)N 


‘11 


1N+1 


dN+ll  '  "  '  dN+lN+l 


1  0 
0  1 


0  0  ...  1 


(4.9) 


If  the  inverse  exists  then  (4.8)  can  be  solved  for  the  0 


m(n) 


m(l) 

W) 


8 


m( N+l) 


11 


N+l  1 


1N+1 


N+l  N+l 


11 


N+l  1 


(4.10) 


As  required  in  (4.6)  the  O  .  .  =  sgn  {0  ,  J  and  this  choice  of  J  ,  . 
n  m(n)  6  m(nj  m(n) 

guarantees  that  zero  belongs  to  the  convex  hull  of  the  set 


{°m(l)Atn(l)  ’  •  ‘  ’ °ra ( N+l ) Am ( N+l )  1  * 


4.2.2  Finding  the  intersection  .no  ,\-.-i  error  ec-.iacions 

Knowing  the  signs  it  now  is  possible  to  find  the  minimax 

approximation  for  the  subsystem  of  N+l  error  equations.  The  only  condi¬ 
tion  that  needs  to  be  satisfied  is  to  determine  where  the  set  of  N+l 


error  equations  intersect.  At  that  point  is  where  the  ninimax  solution 


46. 


is  found.  Since  the  point  being  searched  for  occurs  at  the  intersection 

of  the  N+l  error  equations,  the  magnitudes  of  all  the  error  equations  are 

equal  at  this  point.  If  the  vector  yb  =  (y^,  ....  y^)  is  the  point  desired 

and  if  £* is  the  magnitude  of  each  error  equation  evaluated  at  y\  then 

r  ,  , (y^)  can  be  written 
m(n) 


1  N  1 
am(n)rm(n)^y  ^  bm(n)  ^ ^  m(n)kyk  Cm(n)C 


(4.11) 


m(n)  =  m(l),  ...»  m(n+l) 


where  A  .  ,  =  (a  .  ....... a  .  ,,,).  Rearranging  (4.11)  yields 

m(n)  m(n)l  m(n)N 


°m(n)el  +  J/.Mk  yk  ‘  bm(n)  »(")-®U) , . . .  ,»CS+1)  (4.12) 


and  this  system  of  N+l  equations  in  N+l  unknowns  can  be  represented  in 
matrix  form: 


°m(l)  am(l)N 


am(l)N  C  bm(l) 


(4.13) 


Jm(N+l)  am(N+l)N  .  .  .  3m(N+l)N  yN  bm(N+l) 
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MICROCOPY  RESOLUTION  TEST  CHART 


The  solution  to  y  is  found  by  inverting  the  matrix  on  the  left-hand 
side  of  (4.13). 

_  _  _  — 1 JL  _ 


e1 

yl 

0m(l)  am(l)l  am(l)N 

• 

_ 

•  •  • 

• 

•  •  • 

• 

•  •  • 

yN 

_  — 

°m(N+l)  am(N+l)l  •**  am(N+l)N 

rl 

V-/ 

e 

C11  "•  C1  N+l 

bm(l) 

• 

= 

•  • 

•  • 

• 

bm(N+l) 

CN+1  l***  CN+1N+1 

bm  (N+l) 

- 

“  - 

(4.14) 


where  the  matrix  on  the  right-hand  side  of  (4.14)  is  the  inverse  matrix. 
Each  component  y^  of  the  vector  y^  is  given  by  the  sum. 


y 


n 


N+l 

jJj  CN+1  k  bm(k) 


and  the  magnitude  of  the  error  is 

N+l 

=  I  J2  Clk  bra(k)  I 


(4.15) 


(4.16) 


It  is  important  to  note  that  y  is  not  necessarily  the  global 
minimax  approximation  to  the  system  of  M  equations  in  N  unknowns. 
However,  y^  is  the  minimax  solution  to  the  subset  of  the  N+l  chosen 
equations.  At  this  point  the  algorithm  tests  whether  y*'  is  the  global 


miniraax  solution. 
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Since  the  system  of  linear  equations  has  a  total  of  M  equations 
in  N  unknowns,  it  may  be  determined  whether  y^"  is  a  global  minimax  approx¬ 
imation  by  evaluating  each  of  the  remaining  M-(N+1)  error  equations  at 
the  proposed  solution  y\  If  any  member  of  the  M-(N+1)  equations  has  a 
magnitude  greater  than  £  when  evaluated  at  y^  then  y^  is  not  the  global 
minimax  approximation.  Otherwise,  y^  is  the  solution  and  the  algorithm 
terminates. 


4.2.3  Exchange  Therorem 

In  the  case  where  y^  is  not  the  solution  a  new  solution  vector 
must  be  computed  for  a  different  set  of  N+l  equations.  As  in  the  case  of 
the  ascent  algorithm  for  one  unknown  the  algorithm  will  exchange  one  of 
the  vectors  in  the  original  set  of  N+l  for  the  vector  associated  to  the 
error  equation  which  evaluated  at  y^  has  the  greatest  magnitude.  Calling 
the  new  vector  Am(N+2)  the  problem  is  to  decide  which  of  the  N+l  vectors 
should  be  replaced,  and  how  this  can  be  accomplished  while  simultaneously 
satisfying  the  condition  that  zero  belong  to  the  convex  hull  of  the  new 
set  of  N+l  vectors. 

It  is  known  that  zero  lies  in  the  convex  hull  of  the  original 

set  of  N+l  vectors  {o  ... .  A  .  . ,  ...,  0  .  A  .  Therefore,  by 

m(l)  m(l)  m(N+l)  m(N+l)  J 

(4.7) 


N+l 

-  I 

n=l 


^m(n)^m(n) 


and 


N+l 

l  0  ,  , 

n=l  ra(n) 


1 


(4.17) 


The  condition  that  zero  lie  in  the  convex  hull  of  the  A  .  .  implies 

m(n) 

that  one  of  the  N+l  vectors,  say  A 


can  be  written  as  a  linear 
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combination  of  the  N  other  vectors: 


Am(j) 


N+l 

l  ~ 

n-1 


m(n) 

^m(j) 


m(n) 


(4.18) 


If  it  is  assumed  that  the  N+l  vectors  A  ,  .  span  N  space  then  it 

m(n) 

is  possible  to  express  the  new  vector  A  as  a  linear  combination  of 

m(N+2) 

the  A  ,  .  for  n  =  1,  ...,  N+l  and  some  scalar  coefficients  X  , 
m(n)  m(n) 


Am(N+2) 


N+l 

*m(n)  Am(n) 


(4.19) 


Equation  (4.19)  can  be  rewritten 


N+l 

Am(N+2)  \n(j)Am(j)  +  \n(n)  Am(n) 

n*j 


Xm(j) 


N+l 

l 


n=l 


9m(n) 

9m(j) 


N+l 

n*2_  Xm(n)  Am(n) 

n/j 


N+l 

n=1  ^m(n)  \n(j) 

n^j 


9m(n) 

m(j) 


m(n) 


(4.20) 


so  that  Am(N+2) 
for  n  =  1,  . . . , 


is  now  expressed  as  a  linear  combination  of  the  A 


m(n) 


N+l  excluding  n= j .  Rewriting  (4.20) 
[N+l  0 


°  =  Am(N+2)  + 


V  \  "»(n) 

n=l  »<J>  em(j) 

n^j 


-  X 


m(n) 


®(n) 


(4.21) 


it  is  seen  that  zero  is  expressed  as  a  linear  combination  of  N+l  vectors 


and  that  the  coefficient  of  each  vector  is  positive  if  j  is  chosen  such 
9m(n) 

that  X  q -  X  -  .  > 0,  or  equivalently  X  ..v/9  /.x>X  ,  x/8  ,  . 

m(j)  m(n)  m(j)  m(l)  m(n)'  m(n) 

for  all  n.  This  condition  ensures  that  the  vector  replaced  by  ^ 

is  such  that  zero  remains  in  the  convex  hull  of  the  new  N+l  vectors. 

The  sign  O  associated  with  vector  A  is  simply  the 

m(N+2)  m(N+2) 

sign  of  rm(N+2)  (x)  without  magnitude  brackets  evaluated  at  y  , 


°m(N+2)  =  Sgn  *bm(N+2) 


N 

-  I 

k=l 


m(N+2)k  k 


(A. 22) 


Since  it  is  desired  that  zero  belong  to  the  convex  hull  of  the  vectors 

modified  by  their  signs,  the  modified  coefficients  X  ,  x  and  8  ,  .  will 
J  b  ’  m(n)  m(n) 

now  be  determined. 

To  find  9  ,  first  recall  that  the  matrix  of  the  o  ,  .A  .  x 
m(n)  m(n)  m(n) 

times  its  inverse  equals  the  identity  matrix 


C11  *  *  *  C1N+1 

0m(l)  3m(l)  1  •"  am(l)N 

1  0  ...  0 

#  m 

•  •  • 

•  •  • 

0  1  ...  0 

CN+1  i’“  CN+1  N+l 

°m(N+l)  ara(N+l)  1  "*  3ra(N+l)N 

0  0  ...  1 

L-  — 

- 

- 

(A. 23) 


Thus, 


N+l 


N+l 


^  °m(n)  Cln  “  1  and  ^  Cln  3m(n)k  =  0  k=1 . N 

n»l  7  n=l 


(A.2A) 


which  is  exactly  the  condition  that  zero  belong  to  the  convex  hull  of 


Che  ®mCa)  Am(n)  »h*r*  the  are  the  lifted  Sm(n). 

Now  the  \  ,  .  are  the  coefficients  of  the  vector  A  .  in 
n'n/  m(N+2) 

the  linear  combination 


N+l 

Am(N+2)  Xm(n)Ara(n) 


(4.25) 


which  is  modified  in  matrix  form  as 


Jm(N+2)’  am(N+2)l’  **•’  am(N+2)N) 


am(l)  am(l)l  *•*  am(l)N 


(Xm(l)’,-*‘Xm(N+l)) 


m(N+l)  m(N+l)l  ***  am(N+l)N 


(4.26) 


and  therefore, 


(Xm(l)’,"*Xm(N+l)) 


(°m(N+2),am(N+2)l'””  3m(N+2)N) 


N+l  1 


'1  N+l 


N+l  N+l 


(4.27) 
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Since 


°m(N+2)Am(N+2) 


N+l 

ara(N+2)  ^m(n)Am(n)  ^ 

n=l 


N+l 

n^^°m(N+2)0m(n)^m(n)  ^  ^m(n)  Am(n* 


(4.28) 


the  modified  A  .  .  are  the  a  .„,_.0  .  .A  ,  .  and  thus  the  modified  ratios 
ra(n)  m(N+2)  m(n)  m(n) 

to  be  computed  are  o  ,  .A  .  ./o  .  .C,  =  0  /M  „.A  ,  . /C 

m(N+2)  m(n)  m(n)  m(n)  In  m(N+2)  m(n)  ln 

The  index  which  has  the  greatest  such  ratio  is  the  index  of  the  vector 
which  will  be  replaced. 

Now  a  minimax  approximation  can  be  computed  from  the  new  set  of 
N+l  equations  and  tested  to  determine  whether  it  is  the  global  minimax 
approximation. 


4.3  Manipulation  of  the  Inverse  Matrix 

Two  important  topics  remain.  The  first  topic  is  concerned  with 
the  nature  of  the  inverse  matrix  in  (4.14).  The  second  topic  deals  with 
convergence  of  the  algorithm.  Convergence  of  the  algorithm  will  be 
discussed  in  the  next  section. 

The  algorithm  described  to  this  point  requires  two  matrix 
inversions  for  the  first  iteration  and  one  inversion  for  each  successive 

iteration.  The  first  two  inversions  are  to  determine  the  signs  0  in 

m 

(4.10)  and  to  find  the  intersection  of  the  initial  set  of  N+l  error 

equations  (see  (4.14)).  For  each  successive  iteration  of  the  algorithm 

a  vector  from  the  old  set  of  N+l  vectors  is  exchanged  for  a  new  vector 

A  (see  section  4.2.3).  The  sign  of  the  new  vector  is  established, 

m(N+2; 

and  the  only  operation  to  be  performed  is  the  location  of  the  point  of 
intersection  of  the  set  of  N+l  new  error  equations.  As  in  (4.14)  this 


3 
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can  be  accomplished  through  a  matrix  inversion  of  the  new  set  of  N+l 

vectors  and  their  associated  signs. 

An  algorithm  with  many  matrix  inversions  is  undesirable  and 

should  be  avoided,  if  possible.  The  reasons  for  this  are  because  matf*ces 

can  be  difficult  to  invert  numerically  (due  to  round-off  error,  and  d?ta 

inaccuracies),  matrix  inversion  algorithms  are  time-consuming. 

Since  only  one  row  of  the  matrix  is  exchanged  for  each  iteratfon 

of  the  ascent  algorithm  it  is  possible  to  eliminate  all  of  the  matrix 

inversions  beyond  the  initial  two.  In  general,  if  only  the  Btfl  row  of  a 

square  matrix  (whose  inverse  is  known)  is  altered,  the  changes  in  the 

inverted  matrix  can  be  predicted  without  inverting  the  altered  matrix: 

Assume  that  an  NxN  matrix  has  rows  A  =  (a  ,,..., a  ,)  and  its 

n  nl  nN 

inverse  has  columns  C  =  (C,  ,...,C  ).  Furthermore,  assume  that  the 

n  in  Nn 

matrix  is  altered  by  exchanging  Aa,  the  matrix's  Bth  row,  for  the 

P  < 

* 

N-vector  v.  The  altered  matrix  will  now  have  rows  A  and  the  inverse 

n 

* 

will  have  columns  .  Note  that  the  only  difference  between  the  matrix 

*  th  *  * 

A  and  A  is  that  the  B  row  of  A  is  A^  =  (v^,...,vN). 

It  now  will  be  shown  that  the  Bth  column  of  the  altered  inverre 

c6*t. 

CB*  .  Cg/<v,Cg>  (Z  -25> 


and  that  every  other  column  of  the  altered  inverse  is 


<V’Cn>CB 


(/.  -30) 


n  *  1, . . . ,N+1 
n  i1  8 
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where  <•,•>  denotes  the  inner  product.  It  turns  out  that  the  transforma- 

* 

tion  in  (4.29)  and  (4.30)  is  such  that  the  matrix  C  is  the  inverse  of  the 
*  *  *  *  * 

matrix  A  (i.e.,  A  C  =  C  A  =  I,  for  I  the  identity  matrix) . 


There  are  four  cases  to  take  into  consideration  to  prove  (4.29) 


and  (4.30), 


Proof : 


Case  (i):  for  m  =  n  =  $, 


<A&  ,Cg  >  =  <v,Cg/<v,Ce»  =  <v,Cg>/<v,Cs>  =  1 


Case  (ii):  for  m  =  8  &n  J  0 


<Ae  ,Cn  >  =  <v,Cn-<v,Cn>CB>  -  <v,Cn>  -  <v,CnXv,C6  > 

■  <v,Cn>  "  <v'Cn>[<v’CB>/<v’CB>1  ■  0 


Case  (iii):  for  m  4  8&n  =  g 

<A»*'CB*>  ■  <Am,CB/<V,CB>>  ‘  <An,-CB>/<v-CB> 

■  0/<v,C  >  •  0 


Case  (iv):  for  m  j*  6  &  n  4  8 


<A  *,C  *>  =  < A  ,C  -<v,C  >C  *>  =  <A  ,C  >  -  <v,C  ><A  ,C0*> 
m  n  m  n  nB  m’  n  ’  n  m’  6 


=  <A C  >-<v,C  >  [<A  CR>/<v,C  >] 
m  n  n  m  p  n 

»  <Am’Cn>  "  <v,C  >[0/<v,C>]  =  <A  ,C  >  =  6 
ran  n  n  m  n  n® 


k  k 

Cases  (i)-(iv)  have  shown  that  C  is  the  inverse  of  the  matrix  A  . 

In  the  algorithm,  one  of  the  rows  of  the  matrix  is  replaced 
with  the  vector  v  -  (o^^,  an(N+2)1>  ....  «n(lI+2)lI)  and  the  inner 
products  between  v  and  the  columns  of  the  original  inverse  the 
<vfC  >  need  to  be  computed  to  take  advantage  of  the  above  proof. 


However,  these  inner  products  have  previously  been  computed  in  the  exchange 


portion  of  the  algorithm  and  are  the  X  ,  .  . 

m(n) 


Specifically, 


\n(n)  °m(N+2)^m(N+2)  ’  SnOi)"* 


(4.31) 


(see  (4.27)). 

There  is  one  condition  which  must  be  satisfied  when  manipulating 

jc 

the  inverse  of  the  matrix  in  (4.14).  If  A  is  a  singular  matrix  then  the 
inner  product  <v,  C^>  =  0  and  the  foregoing  calculations  are  not  valid. 

4*4  Convergence 

It  remains  to  be  shown  that  the  ascent  algorithm  converges  to  a 
solution.  Since  there  are  only  a  finite  number  of  subsets  of  N+l  vectors 

r  Mi 

out  of  the  M  (JN+IJ  =  M  !  / (N+l) ! (M-N-l) !  to  be  exact),  the  solution  would 
be  guaranteed  simply  by  testing  all  possible  combinations.  Since  the 
ascent  algorithm  does  not  check  all  possible  combinations,  it  must  be 
shown  that  the  algorithm  does  not  cycle  -  i.e.,  alternate  from  the  current 
solution  to  a  previous  solution.  This  is  shown  by  proving  that  the 
magnitude  of  the  minimax  error  grows  larger  for  each  iteration  of  the 
algorithm. 

Suppose  that  a  certain  step  of  the  algorithm  yields  the  set  of 

vectors  (A^, . . . and  associated  with  this  set  is  a  computed  solution 

vector  y  and  error  e.  Assume  further  that  in  the  next  step  of  the 

algorithm  that  A^  will  be  replaced  by  A^^  t'1us  creating  the  new  set  of 

vectors  {A_, . . . ,A„. _}  .  Associated  with  the  new  set  of  vectors  is  a  new 
2  N+Z 

solution  vector  y'  and  error  E*. 


From  the  choice  of  A^+2  it  is  known  that 


lrN+2<y>l>lr2(y)l  and  lrN+2(y,)l=lr2Cy,)l  • 

Thus  y  j  y*  or  equivalently 

y-yWO. 

For  n  =  2,  N+l  it  is  true  that 

<0nAn>y-y,>  =  [<°nVy>  ~  V  "  l<0nVy,>  "  bn] 

=  °nrn(y)  ~  °nrn(y,)  =  6-6  ’  * 

Thus  by  (4.33)  and  (4.34) 

e  -  e'  +  0 

For  n  =  N+2, 

<0N+2AN+2’  y_y  >  =  ^<C7N+2AN+2’y>  ~  bN+2  ^  <0N+2\+2 

°N+2rN+2^y*  “  °N+2rN+2  ^y ' ^ >E"e ' 

If  c-e’>0  then  it  must  be  so  chat  <a  A  ,y-y’>>0  for  n  =  2,  . 
(4.34)  and  (4.36).  However,  it  will  now  be  shown  that  if  e 
zero  does  not  belong  to  the  convex  hull  of  the  set  {A^,  ...» 


(4.32) 

(4.33) 


(4.34) 


(4.35) 

’y,>  -bN+2] 

(4.36) 

. .,  N+2  by 
-e’>0  then 
An+2)  which 


contradicts  the  condition  that  the  set  have  a  minimax  solution. 
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Equation  (A. 39)  must  be  violated  if  <0  A  ,y-y'»0  for  all  n=2,...,N+2 

n  n 

since  by  (4.37)  the  0^  are  all  greater  than  or  equal  to  zero.  Thus  it 

has  been  demonstrated  that  «J  A  ,y-y'>  cannot  be  greater  than  zero  for 

n  n 

all  9  and  it  raav  be  concluded  that  the  supposition  e-e'>0  is  false, 
n 

Due  to  (4.35)  the  only  remaining  possibility  is  that  e'-£>0  or 
equivalently  £'>£.  Thus,  for  each  iteration  of  the  algorithm  the  value 
of  the  error  must  increase  and  this  fact  proves  the  algorithm's 
effectiveness . 

4.5  Restrictions  on  the  Ascent  Algorithm 

There  is  an  underlying  assumption  which  makes  the  ascent  algorithm 
operational.  It  concerns  the  condition  of  the  set  of  linear  equations, 
the  extent  to  which  these  equations  are  linearly  independent  of  each  other. 
For  the  algorithm  to  work,  it  is  necessary  to  perform  matrix  inversions 
and  to  manipulate  matrix  inverses.  The  algorithm  proceeds  on  the 
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assumption  that  for  each  matrix  and  inverse  matrix  it  works  with  that 
each  one  will  be  non-singular. 

A  non-singular  matrix  is  a  matrix  whose  determinant  is  non-zero 
which  is  equivalent  to  the  condition  that  the  rows  of  the  matrix  be 
linearly  independent  and  that  the  columns  of  the  matrix  be  linearly 
independent.  Therefore,  if  the  algorithm  is  to  operate  on  matrices  and 
their  inverses,  it  must  be  guaranteed  that  the  inverses  exist.  This  is 
equivalent  to  requiring  that  every  subset  of  N+l  from  the  set  of  M 
vectors  be  linearly  independent. 

* 

This  condition  is  known  as  the  Haar  condition  and  the  ascent 
algorithm  will,  in  general,  only  be  successful  if  the  rows  of  the  matrix 
in  (2.2)  satisfy  the  Haar  condition. 

5.  USING  THE  TRANSVERSAL  FILTER  AS  A  CONTINUOUS  TIME  EQUALIZER 
5.1  Statement  of  The  Equalization  Problem 

When  the  system  or  channel  to  be  equalized  has  a  continuous 
time  impulse  response  h(t),  then  the  convolution  between  h(t)  and  the 
impulse  response  of  a  transversal  filter  f(t)  is 

N 

g(t)  =  l  fnh(t-nT)  (5.1) 

n=l 

where  g(t)  is  the  desired  equalized  response,  the  f  are  the  tap  weights 

of  the  filter,  and  T  is  a  unit  delay  (see  Fig.  1.3).  For  h(t)  measured 

unu  gi.c;  specified,  the  objective  of  equalization  is  to  soiv  . 

the  value  of  the  tap  weights,  the  f  .  As  previously  discussed  in 

n 

section  1.3,  (5.1)  has  no  exact  solution  unless  g(t)  can  be  written  as  a 

* 

Cheney,  E.W. ,  Introduction  to  Approximation  Theory,  McGrar  Hill, 


1966. 


linear  combination  of  the  h(t-nT)  for  n=l,...,N.  In  general,  the  solutions 
to  (5.1)  are  approximate  solutions  and  to  gauge  the  quality  of  these  solu¬ 
tions  an  error  function  r(t)  is  established  as  in  (1.8). 

N 

r (t)  =  g(t)  -  l  f  h(t-nT).  (5.2) 

n=l 

Now,  as  discussed  in  section  1.4  the  norm  of  r(t)  will  serve  as  a 
measure  of  closeness  between  the  approximation  £  fnh(t-nT)  and  the  desired 
function  g(t).  The  Chebyshev  norm  will  be  the  specific  distance  measure 
explored  here.  The  Chebyshev  norm  of  r(t)  defined  on  an  interval  [a,b] 
is  defined  as 


N  N 

|jr(t)!|  =  ]|g(t)  -  l  f  h(t-nT)||  =  max  ]g(t)  -  J  f  h(t-nT)  |  (5.3) 

n=l  a<t<b  n=l 


In  Chapter  2  the  discussion  centered  on  how  to  find  the 
Chebyshev  norm  of  an  overspecified  system  of  linear  equations  and  how  to 
minimize  the  norm  of  such  a  system.  This  process  and  its  results  are 
known  as  the  minimax  approximation  and  it  is  desired  to  find  a  minimax 
approximation  to  the  problem  of  (5.1)  by  minimizing  the  norm  of  r(t) 
as  defined  in  (5.3). 

The  path  to  follow  which  leads  to  a  minimax  solution  of  (5.1) 
is  not  readily  apparent  since  the  practical  problem  of  finding  the 
minimum  point  of  ||r(t)||  has  only  been  discussed  in  terms  of  a  system 
of  linear  equations.  Section  5.2  will  present  a  method  for  solving 


(5.1)  in  the  minimax  sense. 
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5.2  Linearization  of  the  Continuous  Error  Function 

The  function  r(t)  can  be  seen  as  an  error  equation  of  a  linear 
equation  in  N  unknowns  if  (5.1)  is  evaluated  at  a  single  point  in  time 


N 

g(t^)  =  l  y^tj-nT)  .  (5.4) 


Equation  (5-4)  is  of  the  form  of  a  linear  equation 


b  = 


N 


l 

n=l 


a  x 
n  n 


(5.5) 


where  the  values  of  b  and  a^  are  determined  by  the  functions  g(t)  and  h(t) 
evaluated  at  the  points  t^  and  t^-nT.  Specifically,  b  =  g(t^)  and 
3n  =  Mtj-nT). 

Any  interval  [a,b]  of  the  time  axis  has  infinitely  many  points 
belonging  to  it  and  each  of  these  points  may  be  evaluated  in  (5.4)  to 
generate  a  linear  equation  in  N  unknowns.  Thus  (5.1)  can  be  considered 
an  overspecified  system  of  linear  equations  where  the  number  of  equations 
is  infinite  and  the  number  of  unknowns  is  N. 

From  section  2.4  it  is  known  that  the  minimax  approximation  is 
the  minimum  point  of  the  function  ||r||  and  that  the  minimum  point  is 
located  by  the  intersection  of  some  subset  of  N+l  error  equations  from 
the  total  set  of  M.  The  object  of  the  next  chapter  will  be  to  locate 
the  N+l  points  on  the  time  axis  which  when  used  in  (5.2)  will  generate 
the  minimax  approximation  to  (5.1). 


6.  REMEZ  ALGORITHM 


6.1  Process  of  the  Remez  Algorithm 

As  demonstrated  in  Chapter  2,  it  is  always  possible  to  find  the 
minimax  approximation  to  a  finite  overspecified  system  of  linear  equations 
in  N  unknowns.  The  minimax  approximation  is  located  by  some  subset  of 
N+l  equations  and  since  there  are  only  a  finite  number  of  subsets  contain¬ 
ing  N+l  equations  the  location  of  the  minimax  approximation  is  guaranteed 
simply  by  testing  all  possible  subsets. 

When  the  set  of  equations  is  infinite,  the  number  of  subsets 
containing  N+l  equations  is  also  infinite.  Thus,  it  is  not  possible  to 
test  all  possible  subsets  of  the  system.  What  is  sought  is  an  algorithm 
which  sequentially  examines  subsets  of  N+l  equations  such  that  the 

sequence  converges  to  the  optimum  subset  of  N+l  equations.  One  algorithm 

* 

which  accomplishes  this  is  called  the  Second  Algorithm  of  Remez  which  is 
analogous  to  the  ascent  algorithm  described  in  Chapter  4. 

In  the  time  domain  the  error  function  r(t)  (as  defined  in  (5.2)) 
is  of  the  form 


N 

r(t)  =  g(t)  -  7  f  h(t-nT)  a<t<b  (6.1) 

n  - 

n=l 

for  some  interval  [a,b].  As  discussed  in  section  5.2,  (6.1)  is  equivalent 
to  an  overspecified  system  of  linear  equations  (with  an  infinite  number 
of  equations)  in  N  unknowns. 

The  Remez  algorithm  begins  by  choosing  an  arbitrary  subset  of 
N+l  points  from  the  time  interval  [a,b].  These  N+l  points  define  a  subset 
of  N+l  linear  equations  in  N  unknowns  (see  section  5.1)  and  it  is  possible 

*Chenev  ,  L.W. ,  Introduction  to  Approximat ion  Theory,  McGraw  Hill, 
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to  find  the  minimax  approximation  to  this  subset  of  N+l  equations  using 
the  ascent  algorithm  as  described  in  Chapter  4. 

The  set  of  tap  weights  {f^,...,f^}  and  the  minimax  error  e  is 

computed  from  the  subset  of  N+l  equations.  For  {f, . f.,}  to  be  the 

X  N 

global  minimax  approximation  to  r(t)  on  [a,b],  the  magnitude  of  r(t)  must 
be  bounded  by  £.  This  condition  may  be  determined  by  inserting  the 
computed  set  {f^,...,f^}  into  the  continuous  expression  for  r(t)  in 
(6.7).  If  |r(t) |>e  for  any  t  belonging  to  [a,b]  then  the  set  {f^,...,f  } 
does  not  constitute  the  global  minimax  approximation  to  r(t). 

The  subset  of  N+l  equations  is  then  replaced  by  a  different  sub¬ 
set  of  N+l  equations.  This  is  accomplished  by  choosing  a  new  set  of  N+l 
points  from  the  interval  [a,b].  The  procedure  differs  here  from  the 
ascent  algorithm  as  described  in  section  4.1  since  more  than  one  equation 
is  replaced  for  each  iteration.  Generally,  all  N+l  equations  are  replaced 
by  N+l  new  equations  in  each  iteration  of  the  algorithm. 

A  graphical  example  of  exchanging  more  than  one  equation  per 
iteration  is  shown  in  Fig.  6.1.  The  ascent  algorithm  described  in 


Figure  6.1  Example  of  Multiple  Exchanges  in  Ascent  Algorithm 
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Chapter  4  determines  the  global  minimax  approximation  by  computing  a  suc¬ 
cession  of  solutions  for  subsets  of  N+l  equations.  Since  Fig.  6.1  is  an 
example  of  a  system  in  one  unknown  the  subsets  will  contain  two  lines  and 
the  point  x^  is  a  relative  minimax  solution  since  it  is  located  at  the 
intersection  of  two  lines  which  have  slopes  of  opposite  sign.  Furthermore, 
when  evaluated  at  x^  there  are  two  lines  of  greater  magnitude  than  the 
magnitude  of  the  lines  which  intersect  at  x^. 

The  Ascent  algorithm  described  in  Chapter  4  proceeds  by  replacing 
one  of  the  two  intersecting  lines  at  x^  by  the  single  line  of  greatest 
magnitude  evaluated  at  x^.  The  intersection  of  the  two  new  lines  locates 
a  relative  minimax  solution  at  X2*  Repeating  this  process  once  again 
locates  the  global  minimax  solution  x^. 

At  x^  in  Fig.  6.1  if  both  intersecting  lines  which  were  used 
in  locating  x^  are  replaced  in  one  iteration  by  the  two  lines  of  greatest 
magnitude  at  x^  the  algorithm  converges  faster.  The  intersection  of  the 
two  new  lines  immediately  locates  the  global  minimax  solution  at  x^  and 
bypasses  the  examination  of  the  relative  minimax  solution  at 

In  the  Remez  algorithm  every  set  of  N+l  equations  is  replaced 
by  a  new  set  of  N+l  equations.  The  algorithm  sequentially  computes  sets 
of  tap  weights  (one  set  for  each  set  of  N+l  equations)  which  in  the 
limit  converge  to  optimum  values.  Because  the  algorithm  must  eventually 
terminate,  the  criterion  for  ending  the  algorithm  depends  on  the  percent 
difference  between  two  consecutive  computations  of  ||r(t)|j.  When  the 
percent  difference  is  smaller  than  some  specified  percentage  the 


algorithm  ends. 
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6.2  Operational  Details  of  the  Remez  Algorithm 

The  framework  of  the  algorithm  has  been  described,  what  remains 
are  the  operational  details.  There  are  three  parts  of  the  Remez  algorithm, 
they  are:  (1)  ariving  at  the  relative  minimax  approximation  for  the  N+l 
equations  generated  by  part  (3);  (2)  evaluation  of  the  merit  of  the  approx¬ 
imation  in  part  (1);  and  (3)  choosing  a  set  of  N+l  points  from  the  time 
interval  [a,b]  for  use  in  part  (1). 

The  first  set  of  N+l  points  can  be  arbitrary,  and  it  is  assumed 
that  the  algorithm  has  a  set  of  N+l  points  t  such  that  a<t. < . . . <t„ , , <b . 
These  points  are  used  to  generate  an  overspecified  system  of  N+l  linear 

equations  in  N  unknowns  by  evaluating  (5.1)  at  each  of  the  points  t  . 

n 

The  minimax  solution  for  f  is  found  by  the  ascent  algorithm 

which  was  described  in  Chapter  4.  The  only  essential  difference  is  that 

the  portion  of  the  ascent  algorithm  which  performs  manipulations  on  the 

inverse  matrix  (see  section  4.3)  is  invalid  since  more  than  one  equation 

is  exchanged  in  each  iteration  of  the  algorithm.  As  in  the  first  step  of 

the  ascent  algorithm  a  set  of  N+l  signs  o  must  be  computed  for  the  N+l 

n 

equations  and  then  the  (N+l)  x  (N+l)  matrix  of  signs  and  vector  coeffi¬ 
cients  is  inverted  to  find  the  rainimax  approximation  to  this  particular 
subset  of  N+l  equations.  An  inversion  of  a  matrix  to  compute  the  signs 
On  is  unnecessary  due  to  an  alternation  property  of  the  minimax 

approximation  for  continuous  time  signals.  This  property  is  described 

* 

by  a  theorem  called  the  Alternation  Theorem  and  the  property  is  a 
consequence  of  the  requirement  that  zero  belong  to  the  convex  hull  of  each 
subset  of  N+l  equations.  Specifically,  the  alternation  theorem  states 
that  the  minimax  approximation  generated  by  the  set  of  points 

* 
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a<_t1<...<tN+1<b  has  the  property  that  r(tn)  =  -r(t  )  for  n=l . N.* 

Knowing  that  the  magnitude  of  the  error  function  will  alternate  in  sign 

at  the  N+l  points  t  is  equivalent  to  knowing  a  priori  that  the  signs 

O  will  alternate  in  sign  since  C  =  sgn  {r(t  )}  . 
n  n  n 

Each  minimax  approximation  computed  by  the  Remez  algorithm  is 
for  a  subset  of  N+l  equations  generated  by  N+l  points  from  the  interval 
[a,b].  For  each  iteration  the  algorithm  tests  how  good  an  approximation 
it  has  computed.  Later  in  this  chapter  it  will  be  proven  that  the  Remez 
algorithm  achieves  linear  convergence.  Using  this  knowledge  it  can  be 
shown  that  each  successive  tap-weight  vector  is  closer  to  the  optimum 
tap-weight  vector  than  the  same  vector  from  the  previous  iteration. 
Termination  of  the  algorithm  can  be  controlled  by  monitoring  the  distance 
between  successive  tap-weight  vectors  or  by  monitoring  the  percent  dif¬ 
ference  between  successive  values  of  ||r(t)||. 

If  the  algorithm  has  finished  computing  the  minimax  approxima¬ 
tion  to  a  particular  subset  of  N+l  equations  and  the  algorithm  does  not 
terminate,  then  it  replaces  the  N+l  points  of  the  past  iteration.  The 
description  of  how  the  algorithm  makes  this  replacement  refers  to 
Fig.  6.2. 

The  central  idea  in  this  part  of  the  algorithm  is  to  exchange 

one  subset  of  N+l  points  for  another  subset  such  that  the  error  function 

evaluated  at  the  new  subset  of  points  has  magnitudes  unbounded  by  the 

minimax  error  of  the  previous  iteration.  It  is  known  that  the  previous 

iteration  generated  a  minimax  approximation  for  the  N+l  points 

a<t,< . . .<t„, , <b.  From  the  characteristics  of  the  minimax  approximation 
—  1  N+l— 

* 

See  section  A. 4. 
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The  alternation  theorem  requires  r(tn)  =  -r(t  j.)  for  n“l,...,N. 

Thus,  for  each  pair  of  indices  (n,  n+1) ,  r(t)  must  have  a  root  in  the 

interval  [ t^ , tQ+^ ] .  For  every  interval  [t  ,t  the  algorithm  then 

locates  this  root  and  labels  it  z  ...  The  set  of  roots  z  also  contains 

n+1  n 

the  endpoints  a  and  b  so  that,  as  can  be  seen  in  Fig.  6.2  abcissa  number  2, 

the  N+2  roots  are  ordered  a  =  z^<...<z^+2  =  b.  For  each  pair  of  roots  the 

interval  [z  ,  z  . ,  ]  contains  a  single  member  of  the  set  {t, , . . . ,t„. „ 
n  n+l  1  N+1 

Using  the  set  of  roots  {z^, —  »z'u+2^  a  test  set  new  Points  will 

be  generated  for  the  next  iteration  of  the  algorithm.  It  is  a  test  set 

because  the  choices  for  the  set  may  be  modified  before  being  passed  to  the 

next  iteration  of  the  algorithm.  The  test  set  is  found  by  choosing  a 

point  T  from  each  interval  [z  ,  z  for  n  =  1,  ....  N+1  such  that 

n  n  n+1 

O  r(T  )  is  a  maximum  in  the  interval.  The  coefficient  is  the  sign  of  the 
n  n 

error  function  evaluated  at  t  (i.e.,  0  =  sign  {r(t  )}  ). 

n  n  n 

Abcissa  number  3  in  Fig.  6.2  shows  the  test  set  {T^, . . . ,TN+^}. 

This  set  is  modified  only  if  the  value  t'  for  which  |r(t')|  is  a  maximum 

in  the  interval  [a,b]  is  not  a  member  of  the  set  {T^, . . . >TN+1) •  This 

condition  is  equivalent  to  ||r(t  )jj>jr(Tn)|  for  all  n.  If  t'  does  not 

belong  to  the  set  (T^, . . . it  is  desired  to  insert  t'  into  the  set 

and  to  remove  an  appropriate  T  .  This  is  done  by  first  inserting  t'  into 

n 

its  appropriate  ordered  position  within  the  set  {T^, . . . ,TN+^} .  Then  one 
of  the  T^  will  be  adjacent  to  t'  and  will  be  such  that  sgn  {r(T^)}  = 
sgn  (r(t')}  .  If  t'  is  less  than  T^  or  greater  than  then  both  T^ 

and  T.,,,  are  considered  adjacent  to  t'.  It  is  the  T  which  satisfies 
these  conditions  which  is  removed  from  the  set.  The  remaining  and  the 
point  t'  constitute  a  new  set  of  t  (see  Fig.  6.2  abcissa  number  4). 
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Figure  6.2  shows  how  the  new  set  of  t  replaces  the  old  set.  Note 
that  Tg  is  the  adjacent  point  to  t'  which  has  the  property  sign  r(T^)  = 
sgn  r(t')  .  The  algorithm  now  can  begin  the  next  iteration  by  finding  the 
minimax  solution  to  the  set  of  N+l  equations  generated  by  the  new  set  of 
N+l  points.  Figure  6.3  shows  a  flow  chart  for  the  Remez  algorithm. 

6.3  Convergence  of  the  Remez  Algorithm 

It  remains  to  be  shown  that  the  Remez  algorithm  converges  to  a 
best  tap-weight  vector  f  *.  A  best  tap-weight  vector  in  the  minimax  sense 
is  one  for  which  the  associated  minimax  error  e*  bounds  the  magnitude  of 
the  error  function  r(t)  when  r(t)  is  evaluated  using  fn*> 

N 

r (t)  =  g(t)  -  l  f  *h(t-nT).  (6.2) 

n=l  n 

For  each  iteration  of  the  algorithm  there  is  an  ordered  set  of 

points 

a£t1< . . .  < tN+1<.b  (6.3) 

for  which  a  minimax  approximation  f  and  associated  minimax  error  e 
has  been  computed.  Thus, 

e  =  | r (tL) |  =  |r(t2)|  =  ...  =  |r(tN+1)|.  (6.4) 

Furthermore,  for  the  same  iteration  a  new  set  of  ordered  points 

act, '<. . .<  t'  <b  (6.5) 

—  1  N+l— 

is  chosen  such  that  for  each  point  t  |r(t  ')!>£  .  Let 

n  n 
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a  =  min  |r(t  ')|  ,  8  =  max|r(t  '  )|  (6.6) 

l<n<N+l  n  l<n<N+l  n 


So  that  by  the  theorem  of  de  la  Vallee  Poussin  and  by  (6.4)  and  (6.6) 


* 

e  <.  ot  e  <.  8  • 


(6.7) 


The  proof  of  the  fact  that  the  Remez  algorithm  converges  to  a  best  tap- 
weight  vector  f  *  is  made  by  showing  that  both  8  and  a  converge  uniformly 
to  e*  . 

Letting  e' ,  a',  and  8'  denote  the  quantities  defined  in  (6.4)  and 
(6.6)  for  the  next  iteration  of  the  algorithm,  the  minimax  approximation 
to  the  set  of  equations  generated  by  (6.5)  is  such  that 


(-l)V  +  l  f.  h(t  ’  -kT)  =  g(t  ')  (6.8) 

.  ,  k.  n  n 


The  error  e'  can  be  solved  for  using  Cramer's  rule. 


g(t') 

h(tjJ-T)  ... 

h(t|-NT) 

1 

h(tj-T)  ... 

h(t[-NT) 

g(t’) 

h(t^-T)  ... 

h(t^-NT) 

-1 

h(t^-T)  ... 

h(t^-NT) 

• 

• 

• 

4 

• 

• 

• 

• 

• 

* 

• 

• 

* 

• 

• 

• 

• 

• 

• 

g(tN+lJ 

^(^4-1“^)  •  •  • 

h(tN;rNT) 

<-nm 

k(tN+l~'^  • -  * 

(6.9) 
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Expanding  each  determinant  along  its  first  column  and  denoting  the  N+l 
minors  for  n  =  1,  . ..,  N+l,  the  expression  for  e'  becomes 


N+l  N+l 

l  V  /  I  M, 


(6.10) 


If  g(t)  -  £  fnh(t-nT)  then  £'  =  0  and  so  g(t)  can  be  replaced  by  r(t)  =  g(t)  - 

l  f  h(t-nT)  in  (6.10).  Furthermore,  since  r(t  ')  =  -r(t'  ,,)  and  since 
n  n  n+l 

the  minors  all  have  the  same  sign*,  e '  becomes 


N+l  N+l 

e’  ’  I  |M„l|r(tn,)|/  I  |Mn|  . 


(6.11) 


Choosing  0  , 
n 


N+l 

l«J/  I  l\l 


(6.12) 


l  0  =  1  and  0<e  <1 

u  n  n 


(6.13) 


for  n  =  1,  ...,  N+l.  By  (6.13) 


£’  =  l  9nlr^tn'^  I—  n>in]r(tn’) 

l<n<N+l 


(6.14) 


Knowing  that  e'x*  allows  the  inequality 


N+l  N+l 

a '  -a  >e ’ -a  =  y  9  I r ( t  * )  I -a  7  6 

^  n‘  n  1  n 

n=l  n=l 


N+l 

=  l  9  (!r(t  *)i-a) 
n=l  n  n 


See  section  A. 4. 


(6.15) 


72.  . 


because  a  '  and  £8^  =  1  by  (6.7)  and  (6.13).  Assuming  that  there  exists 

a  8  such  that  (1-0H  min(0  )  then  by  (6.6) 

l£n<_N$l 

a*  -  a>(l-0)(0-a)>(l-8)(e*-a).  (6.16) 

Since  e*<8(see  (6.7)).  Furthermore,  by  (6.16) 

(e*-a*)  =  (e*-a)-(a’-a)<(e*-a)  -  (1-8)  e*-a)  =  0(e*-a).  (6.17) 


Since  it  is  known  that  O<0<1,  (6.16),  gives  a  relationship  for  the  distance 

from  £*  to  a  for  successive  iterations.  If  the  sequence  of  a’s  are  labelled 

.(1)  _(2) 


,  a  ,  ...  then 


e*  -  a(k)  <0k(e*-a(1))  =  B0k 


(6.18) 


For  B  =  (e*-a^)  a  constant,  the  sequence  a^k^  converges  uniformly  and 
monotonically  to  c*  from  below. 

Furthermore,  since  a'  -  aj>(l-8) (B-a)  by  (6.16)  for  each  term  of 
the  sequence  ,  8^\  ... 


„(k)  *  „(k)  (k)  ,,  Q.-l  r  (k+1)  (k). 

8  -  e  <_  6  -  a  ^  (1-8)  [a  -  a  ] 


<(l-e)_1(6:*-a(k))  <  (1-0)"1  B0k  =  C8k  .  (6.19) 


A 

The  sequence  8  also  converges  uniformly  and  monotonically  to  C  , 
hoever,  it  does  so  from  above.  Thus,  the  algorithm  converges  to  the  best 
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Numerical  Examples  and  Suggested  Future  Research 


7.1  Examples 

Two  examples  are  presented.  The  first  example  shows  a  Gaussian 
function  being  equalized  to  sin(t)/t.  Specifically, 


h(t)  =  .337  exp(-t2/27.6) 

g(t)  =  sin(t) /t.  (7.1) 


Figure  7.1(a)  shows  a  plot  of  h(t)  and  g(t)  in  the  interval  [-3ir,3Tr]. 
Figures  7.1(b)  and  7.1(c)  show  plots  of  g(t)  and  the  equalized  versions  of 
h(t) ,  f&*h  for  a  six  tap  filter  and  f^*h  for  an  eight  tap  filter, 
respectively.  The  tap  spacings  for  £.  and  fn  are  T=it  and  T=3ir/4, 

DO 

respectively. 

The  set  of  tap-weights  computed  for  f,  and  f„  are  shown  with  their 

b  o 

associated  minimax  error  in  Table  7.1. 


e 

.208 

f  (1) 

3.038 

f<2) 

-5.828 

f(3) 

3.783 

f(4) 

3.776 

f(5) 

-5.823 

f  (6) 

3.036 

e 

.021 

f  (1) 

-13.637 

f  (2) 

39.696 

f  (3) 

-52.256 

f  (4) 

25.830 

f<5) 

25.811 

f  (6) 

-52.241 

f  (7) 

39.688 

f  (8) 

-13.635 

TABLE  7.1  Tap-Weights  and  Error  for  f,  and  f_ 
_  b  B 


Note  that  the  error  for  the  eight  tap  filter  is  an  order  of  magnitude 
smaller  than  the  error  for  the  six  tap  filter.  The  number  of  iterations 
to  find  f,  and  f.  was  five. 

D  O 

Figure  7.1(d)  shows  the  error  functions  r^  and  rg  for  both  the  six 
and  eight  tap  filters  when  r(t)  is  evaluated  at  the  values  of  f^  and  fg 
shown  in  table  7.1. 

The  second  example  shows  a  raised  cosine  function  being  equalized 
to  a  triangular  pulse.  In  this  case 

h(t)  =  cos^(t/2) 

g(t)  -  2-4 1 1 1  /it.  (7.2) 

Figure  7.2(a)  shows  a  plot  of  h(t)  and  g(t)  in  the  interval  [-tt,tt]. 

Figures  7.2(b)  and  7.2(c)  show  plots  of  g(t)  and  the  equalized  versions 
of  h(t),  fg*h  for  eight  tap  filter  and  f^**1  for  a  twelve  tap  filter, 
respectively.  The  tap  spacings  for  fg  and  f  ^  are  T  =  tt/4  and  T  =  n/6, 
respectively. 

The  set  of  tap  weights  computed  for  fQ  and  f  are  shown  in 

o  Lc 

Table  7.2  along  with  their  associated  minimax  errors.  Again,  note  that 
the  error  is  substantially  reduced  for  the  longer  length  filter. 


.2' 


e 

f(l) 

.263  £ 

4.281  f(l) 

.200 

7.303 

f(2) 

-12.982 

f(2) 

-16.846 

f(3) 

15.664 

f  (3) 

15.261 

f(4) 

-6.358 

f  (4) 

-14.068 

f(5) 

-6.362 

f(5) 

18.958 

f  (6) 

15.666 

f  (6) 

-9.907 

f<7) 

-12.982 

f<7) 

-9.914 

f(8) 

4.281 

f(8) 

18.962 

f(9) 

-14.069 

f(10) 

15.262 

f(ll) 

-16.846 

f  (12) 

7.303 

TABLE  7.2  Tap-Weights  and  Frror  for  fQ  and  f 

_ _  O _ 1|C 

Figure  7.2(d)  shows  the  error  functions  r0  and  r. 0  for  the 

o  12 

eight  and  twleve  tap  filters  when  r(t)  is  evaluated  at  the  values  of 
f  and  f  shown  in  Table  7.2.  The  number  of  iterations  needed  to 

O  14 

find  fQ  was  four  and  the  number  needed  for  f . .  was  five. 

O  14 


7.2  Future  Research 


The  algorithms  discussed  in  this  work  are  dependent  upon  the 
Haar  condition  being  satisfied  (see  section  4.5).  This  is  not  to  imply 
that  the  minimax  solution  does  not  exist  when  the  Haar  condition  is  not 
satisfied,  but  the  algorithm  fails  to  find  the  solution  if  this  is  the 
case.  An  algorithm  which  successfully  operates  independent  of  the  Haar 
condition  is  necessary  before  minimax  equalizers  can  be  universally 
implemented. 

Investigation  of  the  stability  of  the  Ascent  and  Remez  algorithms 
in  the  presence  of  various  types  of  noise  is  a  topic  of  interest.  Noise 
is  a  practical  problem  and  a  comparaitve  study  of  the  effects  of  noise 
on  minimax  equalizers  with  respect  to  other  types  of  equalizers  is 
important . 

Finally,  there  seems  to  be  a  very  strong  potential  for  minimax 
equalizers  in  adaptive  equalization.  Adaptive  equalization  is  the 
technique  of  varying  the  tap-weights  of  a  transversal  filter  to  compen¬ 
sate  for  variation  in  the  communication  channel  or  system  being  equalized. 
Since  the  Ascent  algorithm  only  works  with  N+l  pieces  of  data  at  a  time, 
the  equalization  of  a  system  or  channel  impulse  response  is  determined 
by  a  particular  subset  of  N+l  pieces  of  the  total  input  data.  An  adaptive 
equalizer  in  the  minimax  sense  promises  to  be  very  fast  for  the  following 
reasons:  the  input  of  new  data  never  increases  the  difficulty  of  the 

•  .  ... 

and  new  pieces  of  data  do  not  necessarily  change  the  current  solution. 

Only  if  the  magnitude  of  the  error  curve  associated  with  the  new  piece 
of  data  has  error  unbounded  by  the  current  solution  is  it  necessary  to 
recompute  the  solution. 
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APPENDIX 


A.l  Proof  that  the  Convex  Hull  Contains  Zero 


For  the  overspecified  system  of  linear  equations  as  defined  in 


(2.2), 


11 


*21 


IN 


2N 


^1  ' •  *MN 


xl~ 

V 

X2 

II 

b2 

.  .  4s  ’ 

1 

> 

(A.l) 


each  row  of  the  matrix  can  be  represented  as  an  error  function 


r  (x)  = 
m 


N 

b  -  y  a  x  | 
m  u.  mn  n 1 
n=l 


(A.  2) 


which  can  be  written  in  inner  product  form, 


r  (x)  =  b  -  <  A  ,x  > 
m  1  m  m’ 


(A. 3) 


where  A  =  (a  ,  a  .  ....  a  )  and  x  =  (x.  ,  ....  x,.) . 

It  will  now  be  shown  that  the  minimax  solution,  y,  of  the 
overspecified  system  in  (A.l)  is  such  that  the  set  of  vectors 

{Am  I  rm(y)  =  ! !r(x) I 


(A.  U) 


contains  zero  in  its  convex  hull. 
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The  norm  of  r(x)  is  defined  as 


'<*>11  -  l<m<M  <'„<*» 


(A. 5) 


and  an  example  of  the  norm  of  r(x)  for  the  case  of  one  unknown  is  shown  as 
the  crosshat ched  portion  of  Fig.  A.l. 


Figure  A.l  Proof  of  the  Convex  Hull  Containing  Zero 


Assuming  y  is  not  the  minimax  solution  it  will  be  shown  that  zero  is  not 
contained  in  the  convex  hull  of  the  set  in  (A. 4)  which  will  constitute  a 
proof  by  asserting  the  contrapositive. 

If  y  is  not  the  solution  (refer  to  Fig.  A.l)  then  y  is  not  the 
minimum  value  of  ||r(x)||  and  for  some  vector  h. 


I  |r(y+h) | |  <  |  |r(y) ! | 


(A. 6) 
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Therefore,  it  has  been  shown  by  contradiction  that  if  the  set  in 
(A. 4)  contains  zero  in  its  convex  hull,  then  y  must  be  the  minimum  point 
of  the  function  | |r(x) j  j . 

A. 2  Property  of  the  Convex  Hull  of  a  Set 

A  set  S  is  convex  if  for  any  two  vectors  y  and  z  which  belong  to 
S  the  linear  combination 


6y  +  (l-8)z  ,  0<6<1 


(A.  14) 


also  belongs  to  S.  It  will  be  shown  that  for  any  set  of  vectors 
{f^,  ...,  f^}  which  belong  to  a  convex  set  S,  the  linear  combination 


N 


l 

n=l 


3  f 


n  n 


for 


N 


l 

n=l 


0 

n 


1,  O<0  <1 

—  n— 


(A. 15) 


also  belongs  to  S. 

Rewriting 


(A. 15) 


n=l 


0  f 
n  n 


N-l 

l  0  f  +  0  f 

n  n  N  N 

n=l 


(A. 16) 


and  defining  L, 


N-l 

l  =  y  e 

n=l  n 


equation  (A. 16)  can  be  rewritten 


N  N-l  8  f 

l  9  f  =  L  l  +  9  f 

n=l  n  n  n=l  L  N  N 


(A. 17) 


(A. 18) 


Since 


9  N-l  6 

0<-^<l,  I  ~j7  =1  »  L  +  0  =1  (A.  19) 

n=l 

then  by  induction  on  N-l,  (A. 18),  is  proved  to  be  a  member  of  the  set  S. 

A.  3  Proof  that  the  Minimax  Solutions  are  Determined  by  Subsets  of  N+l 
Elements 

It  will  now  be  shown  that  the  number  of  elements  contained  in  the 
set  (A. 4)  need  not  be  greater  than  N+l.  From  the  discussion  in  section  A.l 
it  is  known  that  zero  is  contained  in  the  convex  hull  of  the  set, 

{A  |r  (y)  =  | |r(x) | |}  .  (A. 20) 

m  m 

The  proof  will  be  made  by  showing  that  a  vector  Z  which  belongs  to  the 
convex  hull  of  a  set  S  is  expressible  as  a  convex  linear  combination  of 
N+l  or  fewer  elements  of  S. 

Thus,  Z  is  expressed  as 


Z  = 


K 


l 

n=  0 


0  V 

n  n 


(A.  21) 


where 


n 


n=G 


n 


(A. 22) 


and  the  V  belong  to  S.  If  K  is  minimal  then  each  0  in  (A. 22)  is  strictly 
n  n 

greater  than  zero. 


Defining  a  set  of  vectors  Z^, 


Z  =  V  -  Z 
n  n 


(A. 23) 


for  n  =  1,  ...»  K  the  Zr  are  linearly  dependent  because 


K  K  K  K 

l  ez  =  l  8  (v  -z)  =  l  ev  -  z  7  e  =  z-z  =  o 

Lt  «  n  “  T>  n  w  n  m  L  n 


“  n  n  **  n  n 
n=o  n=o 


(A. 24) 


Furthermore,  if  K>N  then  the  Z^  are  linearly  dependent  since  there  are  more 

than  N  vectors  in  a  space  of  dimension  N. 

If  Zq  is  removed  from  (A. 24)  then  suppose  a  set  of  not  all  zero 

coefficients  a  exists  such  that 
n 


y  a  Z  =  0  . 
*■,  n  n 
n=l 


(A. 25) 


Defining  a  =  0,  for  all  X 
o 


K  K  K 

y  (0  +Xa  )Z  =  y  0  Z  +  X  y  a  Z  =  0  +  (0)X  =  0. 
„  n  n  n  n  n  n  n 

n=0  n=0  n=0 


(A. 26) 


Now  X  is  chosen  so  that  |X|  is  a  minimum  and  so  that  one  of  the  coefficients 
in  (A. 26)  disappears,  i.e. 


9  +  Xa  =0. 
n  n 


(A. 27) 


The  remaining  coefficients  in  (A. 26)  are  all  non-negative  since  |X|  is  a 


minimum  and  at  least  one  is  non-zero  since  0  +Xa  =  9  >  0.  Using  this  X 

o  o  o 

a  term  in  (A. 26)  drops  out.  Replacing  Z  with  (A. 23)  in  (A. 26)  yields 

n 


Equation  (A. 28)  implies 


K 

2  l 


n=  0 


K 

(6  +\a  )  =  7(6  +Xa  )V 
n  n  L  n  n  n 
n=  0 


n 


(A. 29) 


which  may  be  written 


K  K 

Z  =  [  7  (0  +Xct  )V  ]  /  [  7(0  +Xct  )]  (A. 30) 

1  L  „  n  n  n  L '  n  n 

n=  0  n=0 

which  is  valid  since  the  coefficient  of  Z  on  the  left-hand  side  of  (A. 29) 
is  non-zero. 

Since  one  of  the  terms  in  (A. 30)  is  zero  Z  has  been  expressed  as  a 
sum  of  K  terms  instead  of  K+l  terms  as  in  (A. 21).  Thus  the  assumed  minimality 
of  K  has  been  shown  false  and  it  must  be  that  K<N.  For  K<N  the  number  of 
terms  in  (A. 21)  is  N+l  or  fewer. 

A. 4  Alternation  Theorem 

For  the  continuous  equalization  case  the  error  function  is  defined 
as  in  (5.2) 

N 

r (t)  =  g(t)  -  l  f  h(t-nT) .  (A. 31) 

n-I  n 

For  the  norm  of  r(t)  to  be  minimized  it  will  be  shown  that  for  a  certain 
set  of  N+l  points  t  <...<t  .  from  the  interval  [a,b]  that 
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For  an  interval  [a,b]  a  set  of  N+l  distinct  points  t  is  a 

1  N+l 

vector  t 

t  —  •••»  (A. 33) 


which  has  an  associated  vector  t 

tn  =  (h(tn-T) . h(tn-Nt& 


(A. 34) 


where  t  in  (A. 34)  is  a  component  of  the  vector  t  in  (A. 33). 

The  Haar  condition  (see  section  4.5)  states  that  every  set  of  N 
vectors  of  the  form  (A.  34)  must  be  linearly  independent.  This  is  the  same 
as  requiring  that  the  determinant 


D[tr 


,t:N+l] 


h(t1-T)  ...  h(trNT) 

h(t2-T)  h(t2-NT) 


h(tN+rT) 


h^tN+l-NT^ 


(A. 35) 


be  non-zero. 

Since  determinants  are  continuous  functions  (determinants  are 
defined  by  polynomial  expansions)  it  can  be  shown  that  all  determinants 

ui  Li' ^ i. ..  s...  -  ^  l ns t 

Dttj*  <  0  <  D  fu2>  •••»  (A. 36) 

for  u,<...<u^+2*  Since  the  determinant  is  a  continuous  operator  there 


exists  a  X  such  that  0<X<1  and 


» *  •  •  > (T— ^^UN+1^  ~  ^  (A. 37) 

If  the  determinant  in  (A. 37)  is  zero  then  for  some  distinct  m  and  n 

Xt  +  (l-X)u  =  Xt  +  (l-X)u  (A. 38) 

n  n  m  m 

or 

X(t  -t  )  =  (1-X) (u  -u  )  .  (A. 39) 

n  m  m  n 


From  (A. 39)  t  -t  is  of  opposite  sign  from  u  -u  which  is  a  contradiction 
n  m  n  m 

of  the  assumption  that  the  t^  and  u^  are  similarly  ordered. 

For  | |r(t) ||  to  be  a  minimum  it  is  known  from  sections  A.l  and 
A. 3  that  zero  must  belong  to  the  convex  hull  of  some  subset  of  N+l  or  fewer 
vectors  {o^t^, . . . ,aN+1tN+1}  •  Thus,  from  section  A. 2 

N+l 


I 

n=l 


0  0  t  =  0 
n  n  n 


(A. 40) 


for 


O<0  <1 
—  n— 


N+l 


l  e 

n=l 


n 


1. 


(A. 41) 


Rewriting  (A. 40) 


/N 


N+l  0  O 

I  11  n 

n=2 


■*,0 


11 


t 

n 


(A. 42) 


and  solving  (A. 42)  using  Cramer's  rule 


92 


Using  the  property  of  determinants  which  states  that  a  column  interchange 
changes  the  sign  of  the  determinant  (A. 43)  becomes 


®L°L 


D[tl»*”’tn-l,tn+l 


»*** 


D[t2»  •  •  • » ] 


(-D 


n-1 


(A. 44) 


and  since  both  determinants  in  (A. 44)  are  similarly  ordered  their  signs 
are  the  same.  Thus, 


sgn 


-8  a  , 

(  5— 1)  -  (-I)*1 
9iai 


and  because  each  0  is  positive, 
n 


(A. 45) 


sgn  an 


(-l)n  sgn  . 


(A. 46) 


Equation  (A. 46)  shows  that  the  must  alternate  in  sign  thus  proving 
(A. 32). 
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A. 5  Program  Listings 

The  following  is  a  listing  of  the  ascent  algorithm  and  the  Remez 
algorithm.  The  subroutines  MINV  and  RIMI  used  in  the  algorithms  are  IBM 
scientific  subroutine  programs  and  are  not  listed  here.  Both  algorithms 
are  programmed  in  Fortran  IV  on  a  DEC- 10  computer. 


A. 5.1  Ascent  Algorithm 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAM  ASCENT 
PURPOSE 

TO  OBTAIN  THE  BEST  SET  OF  TAP-WEIGHTS,  F(J) , 

IN  THE  MINIMAX  SENSE  FOR  AN  OVERSPECIFIED 
SYSTEM  OF  LINEAR  EQUATIONS  GENERATED  BY  CON¬ 
VOLUTION  SUM.  HX.DAT  AND  GX.DAT  ARE  THE  UN- 
EQUALIZED  AND  DESIRED  RESPONSES,  RESPECTIVELY, 
AND  MUST  BE  SPECIFIED  BY  THE  USER. 

USAGE 

.EX  ASCENT,  MINV 

DESCRIPTION  OF  PARAMETERS 

MINV  -SUBROUTINE  WHICH  FINDS  THE  INVERSE  OF 
AN  (N  X  N)  MATRIX  PROVIDED  THE  MATRIX 
IS  NON-SINGULAR. 

REMARKS 

THE  PROCEDURE  ASSUMES  THAT  THE  ROWS  OF  THE 
MATRIX  REPRESENTING  THE  OVERSPECIFIED  SYSTEM 
OF  LINEAR  EQUATIONS  SATISFY  THE  HAAR  CONDITION. 

METHOD 

THE  SOLUTION  IS  FOUND  BY  SUCCESSIVE  SOLUTIONS 
TO  SUBSYSTEMS  OF  N+l  EQUATIONS.  SINCE  THERE 
ARE  A  FINITE  NUMBER  OF  SUCH  SUBSYSTEMS  AND 
SINCE  THE  ERROR  GROWS  FOR  EACH  ITERATION  THE 
ALGORITHM  IS  GUARANTEED  TO  FIND  THE  SOLUTION 
IN  A  FINITE  NUMBER  OF  ITERATIONS. 


INITIALIZATION 


DIMENSION  A(8,2) ,B(3,3) ,THETA(3) ,Y(3) 

DIMENSION  BR(8) ,R(8) , AMBDA(3) ,JS(3) ,JL(3) ,M(3) 
REAL  MU 

DATA  MR, MC, ME/8, 2, 3/ 

CALL  IFILE(21, 'HX.DAT') 

CALL  IFILE(22, 'GX.DAT') 

DO  10  I-l.ME 
JS(I)-I 


10 


C  - 

C  CREATE  SYSTEM  MATRIX 

C  FROM  HX.DAT  DATA  FILE 

C  - 

J=1 

DO  30  1=1, MR 

I F ( I . GT . MR-MC+1 )  GO  TO  20 
READ (2 1,1)A(I,J) 

GO  TO  30 
20  A(I,J)=0. 

30  CONTINUE 

DO  50  1=1, MR 
DO  50  J=  2.,  MC 
K=I-J+1 

IF(K.LE.O)  TO  TO  40 
IF (K.GT. MR-MC+1)  GO  TO  40 
A (I , J)=A(I-1 , J-l) 

GO  TO  50 
40  A(I, J)=0. 

50  CONTINUE 

DO  5? 


C  COMPUTE  SIGNS  WHICH  FORCE  ZERO 

C  TO  LIE  IN  THE  CONVEX  HULL  OF  THE  N+l  VECTORS 

c  - 

SS  READ(2i,j)8Rfi) 

DO  60  1=1,  ME 
60  B(1,I)=1. 

DO  70  K=2 ,ME 
DO  70  L=1,ME 
N=JS (L) 

70  B(K,L)=A(N,K-1) 

IHOLD=ME 

CALL  MINV (B , IHOLD ,DUM , JL , M) 

DO  80  1=1, ME 

80  THETA(I)=B(I,1)/ABS(B(I,1)) 

C  - 

C  INITIAL  COMPUTATION  OF  MINIMAX 

C  SOLUTION  FOR  SUBSYSTEM  OF  N+l  EQUATIONS 

C  - 

DO  90  1=1, ME 
90  B(I , 1)=THETA(I) 

DO  100  K=2 ,ME 
DO  100  L  =1  ,ME 
N=JS(L) 

100  B(L,K)=A(N,K-1) 

IHOLD=ME 

CALL  MINV (B, IHOLD, DUM.JL.M) 

IF(ABS(DUM) .EQ.O.)  WRITE(5,2) 

110  DO  120  K=1 ,ME 

Y(K)=0. 

DO  120  L=1 ,ME 
N=JS(L) 

120  Y(K)=Y(K)+B(K,L)*BR(N) 
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C  - 

C  EVALUATE  M  ERROR  EQUATIONS  AT 

C  REALTIVE  MINIMAX  SOLUTION 

C  - 

DO  140  K=1 ,MR 
R(K)=0. 

DO  130  L=2,ME 

130  R(K)=R(K)+A(K,L-1)*Y(L) 

140  R(K)=R(K)-BR(K) 

C  - 

C  FIND  THE  ERROR  EQUATION  WITH  MAXIMUM 

C  MAGNITUDE  AND  DENOTE  ITS  INDEX  AS  ALPHA 

C  - 

JALPHA=1 

X=ABS(R(1)) 

DO  150  1=2, MR 

IF(X.GT.ABS(R(I))  GO  TO  150 
X=ABS (R(I) ) 

JALPHA=I 
150  CONTINUE 

C  - 

C  HAS  A  GLOBAL  MINIMAX  SOLUTION  BEEN  FOUND? 

C  - 

DO  160  1=1, ME 

160  IF(JALPHA.EQ.JS(I))  GO  TO  220 

IF(X.LE.ABS(Y(1)))  TO  TO  220 

c  - 

C  THE  GLOBAL  MINIMAX  SOLUTION  HAS  NOT  BEEN 

C  FOUND.  EXPRESS  THE  VECTOR  WITH  INDEX  JALPHA  AS 

C  A  LINEAR  COMBINATION  OF  THE  N+l  OLD  VECTORS 

C  - 

MU=R( JALPHA) /X 
DO  180  K=1,ME 
AMBDA(K)  =  0. 

DO  170  L=1,MC 

170  AMBDA(K)=AMBDA(K)+A(JALPHA,L)*B(L+1 ,K) 

180  AMBDA(K)=AMBDA(K)+MU*B(1,K) 

C  - 

C  DETERMINE  WHICH  OF  THE  N+l  OLD  VECTORS 

C  WILL  BE  EXCHANGED  FOR  VECTOR  WITH  INDEX  JALPHA 

C  - 

Q=MU*AMBDA(1)/B(1,1) 

NBETA=1 
DO  190  1=2, ME 

IF(Q ,GT.MU*AMBDA(I) /B(l, I) )  GO  TO  190 
Q=MU*AMBDA(I)/B(1,I) 

NBETA=I 

190  CONTINUE 
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C  - 

C  REVISE  INVERSE  MATRIX 

C  - 

DO  200  K=1,ME 

200  B (K , NBETA) =B (K , NBETA) / AMBDA (NBETA) 

DO  210  J=1,ME 
IF(J.EQ. NBETA)  GO  TO  210 
DO  210  K=1 ,ME 

B(K,J)=B(K,J)-AMBDA(J)*B(K, NBETA) 
210  CONTINUE 

C  - 

C  RE- INITIALIZE 

C  - 

JS (NBETA) =JALPHA 
GO  TO  110 

C  - 

C  INPUT-OUTPUT 

C  - 

210  DO  230  1=1, ME 

230  WRITE(5,3)Y(I) 

1  FORMAT (F) 

2  FORMAT ( IX,' DET  OF  MATRIX  IS  ZERO') 

3  FORMAT ( IX , 1F10 . 8) 

C  ■  - 

C  END  OF  PROGRAM 

C  - 

CALL  EXIT 
END 


A. 5. 2  Reraez  Algorithm 


PROGRAM  REMEZ 
PURPOSE 

TO  OBTAIN  THE  BEST  SET  OF  TAP  WEIGHTS,  F(J) ,  IN  THE 
MINIMAX  SENSE  FOR  THE  APPROXIMATION  G(X)=F(1)H(X-T)+ 
...+F(N)H(X-NT).  THIS  IS  WHERE  H(X)  AND  G(X) 

ARE  DEFINED  FUNCTIONS  AND  THE  COEFFICIENTS 
F(l) , . . . ,F(N)  ARE  TO  BE  DETERMINED.  THE  APPROXIMATION 
IS  OVER  A  USER  DEFINED  INTERVAL  [A,B]  and  T  IS  A  USER 
DEFINED  CONSTANT. 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


USAGE 

-EX  REMEZ , GX , HX , RTMI , MINV 

DESCRIPTION  OF  PARAMETERS 

GX  -SUBROUTINE  WHICH  RETURNS  A  SINGLE  OUTPUT 
VALUE  FOR  A  SINGLE  INPUT  VALUE  ACCORDING 
TO  THE  FUNCTION  G(X).  THE  FIRST  LINE  OF 
SUBROUTINE  GX  MUST  BE: 

SUBROUTINE  GX(Y,Z) 

WHERE  Y  IS  THE  INPUT  VALUE  AND  Z  IS  THE 
RETURN  VALUE. 

HX  -SUBROUTINE  WHICH  RETURNS  A  SINGLE  OUTPUT 
VALUE  FOR  A  SINGLE  INPUT  VALUE  ACCORDING 
TO  THE  FUNCTION  H(X).  THE  FIRST  LINE  OF 
SUBROUTINE  HX  MUST  BE: 

SUBROUTINE  HX(Y,Z) 

WHERE  Y  IS  THE  INPUT  VALUE  AND  Z  IS  THE 
RETURN  VALUE. 

RTMI  -SUBROUTINE  WHICH  WHEN  GIVEN  A  FUNCTION  R(X) 
FINDS  A  ROOT  BETWEEN  TWO  VALUES  OF  THE 
FUNCTION  R(X1)  AND  R(X2) .  R(X1)  AND  R(X2) 

MUST  BE  OF  OPPOSITE  SIGNS. 

MINV  -SUBROUTINE  WHICH  TAKES  A  (N  X  N)  MATRIX  AND 

FINDS  ITS  INVERSE  ASSUMING  THAT  TKE  DETERMINANT 
OF  THE  MATRIX  IS  NONZERO. 

REMARKS 

THE  PROCEDURE  ASSUMES  THAT  THE  GENERATED  FUNCTIONS 
H(X-T) , . . . ,H(X-NT)  SATISFY  THE  HAAR  CONDITION.  THE 
CONDITION  REQUIRES  THAT  EVERY  SET  OF  N  DISTINCT 
VECTORS  H(X(I)-T),...,H(X(I)-NT)  FOR  P1,...,N 
BE  INDEPENDENT. 

SUBROUTINES  USED 

GX  AND  HX  ARE  SUBROUTINES  WHICH  MUST  BE  WRITTEN 
BY  THE  USER.  THEY  ARE  INPUT  SUBROUTINES. 

RTMI  AND  MINV  ARE  STANDARD  SCIENTIFIC  SUBROUTINES 
WHICH  CAN  BE  FOUND  ON  MOST  COMPUTING  SYSTEMS. 


METHOD 

SOLUTION  OF  APPROXIMATION  FOR  MINIMAX  COEFFICIENTS 
IS  ACCOMPLISHED  BY  AN  ITERATIVE  PROCESS  KNOWN 
AS  THE  REMEZ  ALGORITHM.  THE  ALGORITHM  SUCCESSIVELY 
MINIMIZES  THE  LARGEST  N+l  DISCRETE  ERRORS  ON  THE 
INTERVAL  [A,B]  UNTIL  IT  CONVERGES  TO  A  BEST  SOL- 
TION.  FOR  REFERENCE  SEE:  INTRODUCTION  TO 
APPROXIMATION  THEORY,  E.W.  CHENEY,  MCGRAW  HILL  1966. 
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C 

C  INITIALIZE  PROGRAM 

C 

DIMENSION  X(9) ,G(9) ,R(9) ,SIGMA(9) ,Z(10) ,H(9,8) 

COMMON  TAU,N,F(9) ,ME,MC,MR 
REAL  NORM 
EXTERNAL  RX 

DATA  N/8/,MR,MC,ME/9,8,9/ 

DATA  EPS/. 001/, IEND/100/, STOP/O./ 

CALL  0FILE(21, 'GDAT.DAT* ) 

CALL  OFILE(22,' ADAT.DAT') 

CALL  OFILE(23, 'RDAT.DAT') 

CALL  OFILE( 24, ’HDAT.DAT') 

£  INPUT  PARAMETERS,  CREATE  INITIAL  X-VECTOR,  AND  COMPUTE 

£  VALUES  IN  FUNCTION  MATRICES 

C 

WRITE (5,1) 

READ(5 ,2)A,B ,TAU 
SINT-(B-A)/ (N+2) 

DO  10  1=1, N+l 
AI=I*SINT 
10  X(I)=A+AI 

20  WRITE(5 ,3) 

WRITE (5 ,5 )X 
DO  30  1=1, N+l 
XI=X(I) 

CALL  GX(XI,GXI) 

G(I)=GXI 
DO  30  J=1,N 
CNTRJ= (2* J-N-l) / 2 . 

XS=XI-CNTRJ*TAU 
CALL  HX(XS ,HXS) 

30  H(I , J)=HXS 

C 

C  COMPUTE  F-VECTOR  AND  THE  INCURRED  ERRORS 

C  CALL  CHEBY(G,ii,F) 

WRITE(5 ,4) 

WRITE(5 ,5)F 
DO  40  I'=1,N+1 
XI=X(I) 

40  R(I)=RX(XI) 

WRITE (5, 6) 

WRITE(5 ,5)R 
C 

C  HAS  A  SOLUTION  BEEN  FOUND? 

C 

TEST=1.-(ST0P/F(1)) 

IF(ABS(TEST) .LT. . 001)  GO  TO  50 
ST0P=F(1) 

GO  TO  70 

50  TINT=(B-A)/50. 

WPITE(5 ,5) 

WRITE(5 ,5) 


WRITE(5 ,9) 

WRITE (5, 5) 

DO  60  K=l,51 
V=A+(K-1) *TINT 
CALL  GX(V,GV) 

CALL  HX(V,HTV) 

RV=RX(V) 

HV=GV-RV 

WRITE (  5 , 7  )  V ,  GV ,  HV ,  RV 
WRITE(21,7)V,GV 
WRITE(22,7)V,HV 
WRITE(23,7)V,RV 
WRITE(24,7)V,HTV 
60  CONTINUE 

GO  TO  210 

C 

C  LOCATE  ZEROES  BETWEEN  ERROR  VALUES 

C 

70  DO  100  1=1, N 

XI=X(I) 

XIPO=X(I+l) 

HOLD=IEND 
SAVE  =£PS 

80  CALL  RTMI (ZIPO,DUM,RX,XI ,XIPO,EPS , IEND, IER) 
IF (IER-1) 100 , 90 ,100 

90  IEND=2*IEND 

EPS=EPS*2. 

IF (IER. EQ .2 )WRITE(5 ,81) 

81  FORMAT (IX, 'ERROR') 

GO  TO  80 

100  Z(I+l)=ZIPO 

IEND=HOLD 
EPS=SAVE 
Z(1)=A 
Z(N+2)=B 
WRITE(5 ,8) 

WRITE(5 ,5)Z 

C 

C  EDIT  X=VECTOR 

C 

DO  110  1=1, N+l 
XI=X(I) 

U=RX(XI) 

L10  SIGMA(I) =U/AB  S (U) 

NORM=0 . 

DO  130  1=1, N+l 
ZI=Z(I) 

ZIPO=Z(I+l) 

STEP  (ZTPO-ZD/IOOO. 

HOLD=0. 

PMAX=ZI 

DO  130  J=l,1000 


T=ZI+J*STEP 

VAL=RX(T) 

SVAL=VAL*SIGMA(I) 

IF (SVAL.LE.HOLD)  GO  TO  120 

HOLD=SVAL 

PMAX=T 

120  IF ( ABS (VAL) , LE . NORM)  GO  TO  130 

Y=T 

NORM=ABS (VAL) 

130  X(I)=PMAX 

C 

C  RE-EDIT  X-VECTOR 

C 

DO  140  K=1,N+1 
XK=X(K) 

RXVAL=RX(XK) 

140  IF (NORM. LE. ABS (RXVAL))K1=1 

IF(Kl.EQ.l)  GO  TO  20 
X1=X(1) 

IF(Y.GT.Xl)  GO  TO  160 
DO  150  K=2 ,N+1 
J=N+3-K 

150  X(J)=X(J-1) 

X(1)=Y 
GO  TO  20 

160  DO  170  1=1, N+l 

IF(Y. LT. Z(I) )  GO  TO  190 
170  IF (Y.L1  X(I))  GO  TO  200 

DO  180  K=1,N 
180  X(K)=X(K+1) 

X(N+1)=Y 
GO  TO  20 
190  X(I)=Y 

GO  TO  20 
200  X(I-1)=Y 

GO  TO  20 
C 

C  INPUT-OUTPUT 

C 

1  FORMAT (IX, 'PLEASE  INPUT  A , B , TAU ' ) 

2  FORMAT (3F) 

3  FORMAT (IX,1 X-VECTOR: ’) 

4  FORMAT ( IX, 'F- VECTOR: ') 

5  FORMAT ( IX, 10F10. 3) 

6  FORMAT ( IX, 'R(T)  EVALUATED  AT  X:’) 

7  FORMAT ( IX, 3F10. 3 ,6X, 1F10. 3) 

8  FORMAT (IX, 'Z-VECTOR: ') 

9  FORMAT (6X, 'X  G(X)  H(X)*F(X) 

C 

C  END  OF  PROGRAM 

C 

210  CALL  EXIT 

END 


ooo  CM  OUCJOU 


101 


C 


COMPUTE  THE  ERROR  FUNCTION:  R=G-<H,F> 


FUNCTION  RX(XIN) 
COMMON  TAU,N,F(9) 
XOUT=0. 

DO  220  J=1,N 
CNTRJ=(2*J-N-l)/2. 
B  =XIN-CNTRJ*TAU 
CALL  HX(B,C) 

20  XOUT=XOUT+C*F(J+l) 
CALL  GX(XIN,D) 
RX=D-X0UT 
RETURN 
END 


FIND  THE  MINIMAX  SOLUTION  TO  SET  OF  N+l  LINEAR 
EQUATIONS  IN  N  UNKNOWNS 


SUBROUTINE  CHEBY (BR,A,Y) 

DIMENSION  A(9,8),B(9,9),TNETA(9),Y(9),BR(9),JL(9),M(9) 

COMMON  TAU,N,F(9), ME, MC, MR 

M0=0 

230  DO  240  1=1, ME 

240  B(I,1)=(-1)**(I+M0) 

DO  250  K=2 ,ME 
DO  250  L= 1 , ME 
250  B(L,K)=A(L,K-1) 

IHOLD=ME 

CALL  MINV(B, IHOLD.DUM, JL,M) 

DO  260  K=1,ME 
Y (K)=0 . 

DO  260  L= 1 , ME 

260  Y(K)=Y(K)+B(K,L)*BR(L) 

IF(Y(l) .GE.O.)  GO  TO  270 
MO=l 

GO  TO  230 
270  RETURN 

END 


